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COMPUTATION  OF  UNSTEADY  TURBULENT  BOUNDARY  LAYERS  WITH  FLOW  REVERSAL 
AND  EVALUATION  OF  TWO  SEPARATE  TURBULENCE  MODELS 

by 

Tuncer  Cebeci*  and  Lawrence  W.  Carr 


SUMMARY 

Recently  a  new  procedure,  which  solves  the  governing  boundary-layer 
equations  with  Keller's  box  method,  has  been  developed  for  calculating 
unsteady  laminar  flows  with  flow  reversal  [1],  In  regions  where  the  stream- 
wise  velocity  contains  flow  reversal,  the  solution  scheme  was  modified  by  a 
procedure  which  accounted  for  the  downstream  influence.  With  this  modifica¬ 
tion,  the  unsteady  flow  over  a  circular  cylinder  started  impulsively  from 
rest  was  successfully  calculated  to  values  of  time  and  space  greater  than 
in  any  previous  solutions.  An  examination  of  unsteady  separation  for  laminar 
flow  was  made  and  revealed  that  the  unsteady  boundary  layer  for  that  flow, 
even  at  large  times,  was  free  of  singularities. 

In  this  report  we  extend  the  method  of  ref.  [1]  to  turbulent  boundary 
layers  with  flow  reversal.  Using  the  algebraic  eddy  viscosity  formulation  of 
Cebeci  and  Smith  [2],  we  consider  several  test  cases  to  investigate  the 
proposition  that  unsteady  turbulent  boundary  layers  also  remain  free  of 
singularities. 

Since  the  solution  of  turbulent  boundary  layers  requires  a  closure 
assumption  for  the  Reynolds  shear-stress  term  and  the  accuracy  of  the  solutions 
depend  on  this  assumption,  we  also  perform  turbulent  flow  calculations  by 
using  the  turbulence  model  of  Bradshaw,  Ferrlss  and  Atwell  [3];  we  solve  the 
governing  equations  for  both  models  by  using  the  same  numerical  scheme  and 
compare  the  predictions  with  each  other,  restricting  the  comparisons  to  cases 


The  study  reveals  that,  as  In  laminar  flows,  the  unsteady  turbulent 
boundary  layers  are  free  from  singularities  but  there  Is  a  clear  Indication 
of  rapid  thickening  of  the  boundary  layer  with  Increasing  flow  reversal. 

The  study  also  reveals  that  the  predictions  of  both  turbulence  models  are  the 
same  for  all  practical  purposes. 
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I.  INTRODUCTION 


The  prediction  of  unsteady  turbulent  boundary  layers  with  flow  reversal 
is  of  Importance  In  a  number  of  aerodynamic  problems,  notably  in  dynamic  stall, 
buffeting  and  gust  studies.  However,  some  of  the  more  popular  turbulence 
models  implicitly  assume  that  the  wall  shear  is  positive  and  their  extension 
to  unsteady  flows  with  flow  reversal  is  not  easy.  It  requires  modifications 
to  the  functional  form  of  the  law  of  the  wall  and  to  the  manner  In  which  the  wall 
shear  is  determined.  Two  near-wall  assumptions  are  considered  here.  In  the 
first,  the  near  wall  grid  point  is  located  in  the  logarithmic  region  and  the 
law  of  the  wall  is  used  to  link  the  flow  properties  at  this  grid  point  to  the 
wall.  In  the  second,  a  Van  Driest  formulation  due  to  Cebeci  and  Smith  [21  is 
used;  this  impl ies  that  the  grid  point  closest  to  the  wall  will  occur  in  the  vis¬ 
cous  sublayer. 

A  further  aspect  of  these  flows  of  current  interest  is  the  possibility  of 
a  singularity  occurring  in  the  reversed-flow  region.  Examples  of  this  phen¬ 
omenon  have  also  been  reported  in  laminar  flows  but,  in  earlier  studies  Cebeci 
[1,4]  and  Bradshaw  [5]  have  shown  that  the  occurrence  is  not  a  feature  of  the 
governing  equations  but  is  due  to  the  limitations  of  the  numerical  procedure 
used.  We  shall  demonstrate  that,  for  the  examples  we  study, there  is  no  indica¬ 
tion  of  such  a  singularity  in  turbulent  flow  either  but  there  is  a  clear  indi¬ 
cation  of  rapid  thickening  of  the  boundary  layer. 

In  addition  to  the  examination  of  wall  functions,  we  have  also  considered 
two  turbulence  models  for  unsteady  flows  without  flow  reversal.  The  algebraic 
eddy- viscosity  formulation  of  Cebeci  and  Smith  (CS)  is  compared  with  the  trans¬ 
port  model  of  Bradshaw,  Ferriss  and  Atwell  [3]  (BF).  Calculations  were  performed 
to  determine  whether  the  representation  of  unsteady  flows  with  strong  pressure 
gradients  requires  that  account  be  taken  of  transport  of  turbulence  quantities. 

As  will  be  shown,  the  predictions  with  both  models  are  nearly  identical  for 
both  steady  and  unsteady  flows  with  and  without  strong  pressure  gradient. 

The  report  has  been  prepared  with  six  main  sections  describing, 
respectively,  the  governing  equations,  the  numerical  procedure,  the  results, 
concluding  remarks,  references  and  the  computer  program  which  uses  only  the 
CS  model. 
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II.  GOVERNING  EQUATIONS 


The  continuity  and  momentum  equations  can  be  written  for  two-dimensional 
unsteady  Incompressible  laminar  or  turbulent  thin  shear  layers  as: 


3U  +  3V 

3x  3y 


0 


0) 


3U 

3t 
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9ue 

u  — -  + 
e  3x 


3t 

ay 
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Here 


T  =  v  |p--  urvr  (3) 

and  we  recall  that  u'  and  v1  denote  fluctuations  about  the  ensemble- 
average  velocity;  u'  and  v'  are  zero  in  unsteady  laminar  flow,  and 
v  au/ay  is  negligible  outside  the  viscous  sublayer  in  a  turbulent  flow. 

These  equations  are  subject  to  the  usual  boundary  conditions,  which  in  the 
case  of  boundary  layers  are 

y  =  0,  u  =  v  =  0;  y  -*■  6  u  ug(x,t)  (4) 

The  presence  of  the  Reynolds  stress  term,  -u  V  introduces  an  additional 
unknown  to  the  system  given  by  Eqs.  (2)  to  (4).  In  this  report  we  present 
calculations  using  two  different  turbulence  models.  One  is  an  algebraic 
eddy-viscosity  formulation  developed  (for  steady  flows)  by  Cebeci  and  Smith 
and  the  other  is  a  transport-equation  model  developed  by  Bradshaw,  Ferriss 
and  Atwell.  In  the  CS  model,  we  write  Eq.  (3)  as 

T  ■  <v  +  em>  ay  (5) 


with  two  separate  formulas  for  em<  In  the  so-called  inner  region  of  the 
boundary  layer  (em)^  is  defined  by  the  following  formula: 


where 


{0.4y[l  -  exp(-y/A)]}2 
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In  the  outer  region  em  is  defined  by  the  following  formula 


<em>o  =  °*0168  /  (ue  ~u)dy 


The  boundary  between  the  Inner  and  outer  regions  Is  established  by  the  contin¬ 
uity  of  the  eddy- viscosity  formulas. 

In  the  BF  model,  which  Is  used  only  outside  the  viscous  sublayer,  we 
assume  t  =  -i? v'  and  write  a  single  first-order  partial -differential  equa¬ 
tion  for  it;  the  equation  was  originally  developed  from  the  turbulent  energy 
equation  but  can  be  equally  well  regarded  as  an  empirical  closure  of  the  exact 
shear-stress  transport  equation.  This  reads 

3/2 

Bt  =  8t  +  u  +  v  9y‘  2ar  ay Sf  <TV  “2al  T"  (9) 

Here  a^  is  a  dimensionless  quantity,  is  a  velocity  and  L  is  the  dis¬ 
sipation  length  parameter,  specified  algebraically  by  L/6  =  f(n)  with 
n  *  y/6  and  f(n)  given  from  an  analytic  fit  to  an  empirical  curve  by 


f(n)  =  0.095  —  0.055(2n  —  1/ 

0.016  exp[-10(n  -  1.1)] 


n  <  0.18 

0.18  £  n  <  1.1 
n  >  1.1 


In  a  more  advanced  version  of  this  turbulence  model  [6]  L  itself  is  deter¬ 
mined  from  a  transport  equation. 

_ _ _  _ t  _ _ _ 

The  turbulent  transport  velocity  VT,  nominally  (p'u1  +  u'v'  )/u’V 
is  proportional  to  a  velocity  scale  of  the  large  eddies  and  is  chosen  to  be 


V  =  2a 


1  u. 


where  g(n)  is  given  by 


T 


g  (n)  = 


n  <  0.5 
0. 5  <  n  <  1.0 


(12) 


Here  is  evaluated  from  the  continuity  equation  (1),  and  a*  is  evaluated 
from  (1)  and  (2)  on  the  assumption  that  the  velocity  u  is  given  by 


for  0  <  y  <  y-j ;  Eq,  (13)  is,  of  course,  a  special  case  of  (16).-*  The  evalua¬ 
tion  of  a*  is  discussed  in  Ref.  [7];  the  last  term  in  (15)  can  be  as  large 
as  half  the  second  (pressure-gradient)  term.  In  unsteady  flow  without  flow 
reversal,  we  use  the  same  inner  "boundary"  conditions  at  y-j  =  50v/u^,  but 
because  of  the  presence  of  the  time-dependent  term  in  (2),  a*  becomes  more 
complicated.  If  we  again  assume  that  (16)  holds  -  remember  that  the  turbulence 
structure  of  the  inner  layer  is  unlikely  to  be  affected  unless  the  external- 
stream  frequency  is  very  high  -  then  (1)  and  (2)  give 


Integrating  we  can  write 
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T  •  T» +  It  /' udy  -  u(j,i>  vr-  *  k  /'  ^ 


0 


08a) 


because 


f  f  i-  dy+  h  0 
3X  J  U  J 

0  T 


We  can  also  write  (18a)  as 
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T  =  X 


»  O  „  .  .  vy.  3u  3u  '1  „  2.  , 

w+3?yi  +  v3t/  irdy  +  uTF(yi}  ~r #  + v  3T  /  (r)dy  (18b) 

n  *  U_  n  T 


or  as 


3U  '  3t 

T1  =  Tw  +  ylF  3T+  “*yl  3T  +  3^yl 


08C) 


where  F  *  u/ut  at  y  =  y-|  and  a*  conies  from  the  last  term  in  (18b)  and 
is  the  same  as  in  steady  flow. 


Equation  (18c)  now  replaces  (15). 


III.  SOLUTION  PROCEDURE 


We  use  Keller's  two-point  finite-difference  method  (called  the  Box 
method)  to  solve  the  system  of  equations  described  in  the  previous  section. 
The  application  of  this  method  to  unsteady  flows  with  no  flow  reversal  using 
the  CS  model  has  been  described  in  Ref.  [8].  Its  application  to  steady  two- 
dimensional  flows  using  the  BF  model  is  described  in  Ref.  [7].  Here  we 
present  a  description  of  the  extension  of  the  CS  model  to  unsteady  two- 
dimensional  turbulent  flows  with  flow  reversal  as  well  as  a  description  of 
the  extension  of  the  BF  model  to  unsteady  turbulent  flows  with  no  flow 
reversal . 


3.1  CS  Method  with  and  without  Flow  Reversal 
As  in  previous  studies  (see,  for  example  [8]),  we  transform  the  equations 

with 

x  =  x/L,  t  =  tufl/L ,  n  =  (u^vx^y  (19a) 

and  a  dimensionless  stream  function  f(>T,n.t’)»  where 

ip  =  (u0vx)^f(x,n,t)  (19b) 

Here  uQ  is  a  reference  velocity,  L  a  reference  length,  and  ip  is  the 
usual  definition  of  the  stream  function  corresponding  to  the  continuity  equati 
(1).  With  the  relations  defined  by  (19)  and  with  the  definition  of  eddy 
viscosity,  equations  (1)  to  (3)  and  the  boundary  conditions  can  be  written 
as 

(bf")'  +  j  ff"  +  m-  =  x  If'  ~~  -  f "  (20) 

6  J  \  ax  3x  3t  / 

n  =  0,  f  =  f'  =  0;  n  +  »  f'  =  ue/uQ  h  (21) 

Primes  denote  differentiation  with  respect  to  n  anJ 


For  simplicity,  we  shall  now  drop  the  bars  on  x  and  t. 
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We  use  two  separate  solution  procedures  to  solve  the  system  given  by 
Eqs.  (20)  and  (21).  When  there  is  no  flow  rev^sa'i  across  the  layer,  we 
us*'  the  standard  Box.  On  the  othp,  liand,  when  there  i_s  flow  reversal,  then 
we  use  the  so-called  zig-zag  Box  as  described  below. 


To  solve  Eqs.  (20)  and  (21)  by  the  standard  Box  method,  we  first  write 
Eq.  (20)  in  terms  of  three  first-order  equations  by  introducing  new  depend¬ 
ent  variables  u(x,n,t),  v(x,n,t),  that  is, 

f1  =  u  (23a) 

(23b) 


u1  =  v 


o>*>-  ♦ 


ax  at' 


(23c) 


We  next  consider  the  net  cube  shown  in  Fig.  1  and  write  difference  approxi¬ 
mations  to  Eqs.  (23).  Equations  (23a, b)  are  approximated  using  centered  dif¬ 
ference  quotients  and  averaged  about  the  midpoint  (x. ,  t  ,  n.  ,  ).  The 

•  n  j 

difference  quotients  which  are  to  approximate  (23c)  are  written  about  the 

midpoint  (x.  ,  ,  t  ,  ,  n-  ,  )  of  the  cube  whose  mesh  widths  are  r . ,  k  , 
l-*5  n—2  J“^  in 

and  h..  This  procedure  yields  the  following  equations: 

J 

(24a) 
(24b) 

(b>)^-v>^ ♦ 1  <*>&  -  ^  *  vft *  vj:"> 

-  <24c> 


■i  ,n 
j 

fi  ,n 
J-l 

-  vft  ■ 

i  ,n 

i  ,n 

-  vft  ■ 

'j 

"Uj-1 

where 


.  .  -ii- 

i  x.  —  x. 


en  = 


xi-^ 


tn"1 

2.(4)  ,  ai 

=  «J(u  )j X  ~  2(u  ““ 


m. 


•ft 


”4  *  f\-l  ~  2fi-l 


1-1 1  -TV*  +  Z6A2-h  -  2Vi>  -  hj'1t(bv,j234 


-(6v)f4i-}(f.)224-4(P3),;;; 
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Fig.  1  Net  cube  for  the  standard  Box  method. 


Here  by  v?34  we  mean  v]"1  ,n  +  +  v]’n_1,  the  sum  of  the  values 

of  Vj  at  three  of  the  four  comers  of  the  face  of  the  box.  Also 

Vi  -  7  ^j"1  +  f3n:i’i_1) 


un-l  2  luj ' 


*  7  <t"  +  f]-?» 


The  resulting  algebraic  system  given  by  Eqs.  (24)  together  with  the 
boundary  conditions,  which  now  become 

fo  *  %  =  O'  uo  s  “e  (25> 

are  nonlinear.  We  use  Newton's  method  to  linearize  the  system  and  solve  the 
linear  system  by  the  block  elimination  method  discussed  in  ref.  [9]. 
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When  there  Is  flow  reversal  across  the  boundary  layer  at  some  x  and  t, 

we  modify  the  standard  Box  method  used  for  Eq.  (23c)  but  retain  that  for 

(23a, b)  and  still  center  them  at  (x^ ,  tn,  n^).  To  write  the  difference 

approximations  for  the  Box  centered  at  (x^  ,,  tn_^,  nj.^)  we  examine 

previously  computed  values  of  ul *J7-  If  up?  >  0,  then  we  use  the  standard 
in  ” 

Box  method:  If  u. V  <  0,  then  we  write  (23c)  for  the  Box  centered  at  P 

J 

(see  Fig.  2)  using  quantities  centered  at  P,  Q,  and  R,  where 


p  "  (x1’  Vv 


Q  E  ^i-V  V  nW 
) 


R  =  (xi+v  Vr  nj-v 

Equation  (23c)  can  then  be  written  as 

(bv)'  (P)  +  \  (fv)(P)  =  x(P)  (eu(Q)  (Q)  +  <pu (R)  (R)  -  ev(Q)  (Q) 


3f 


—  d>v(R)  f£-(R)  +  |^-  (P)D 


Here 


e  = 


xi+l  "  xi 


M+l 


-  x 


i-1 


4>  = 


X1  ~  xi-l 
xi+l  "  x1-l 


(26) 


(27) 


(28) 


Fig.  2  Finite  difference  molecule  for  the  Zig-Zag  differencing. 
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The  resulting  algebraic  system  Is  again  nonlinear  and  its  solution  is  obtained 
by  using  the  procedure  followed  In  the  standard  Box  method. 

3.2  BF  Method  with  no  Flow  Reversal 


The  solution  of  the  governing  equations  for  unsteacty  flows  with  the  BF 
model,  even  with  no  flow  reversal  across  the  boundary  layer,  is  much  more  dif¬ 
ficult  than  with  the  CS  model.  This  is  because  of  the  hyperbolic  nature  of 
the  governing  equations,  together  with  the  nonlinear  boundary  conditions,  which 
play  an  important  role  in  the  solution  procedure.  As  is  common  in  most  (if  not 
all)  methods  that  use  boundary  conditions  away  from  the  "wall,"  the  wall  shear 
stress  is  also  an  unknown  parameter;  it  can  be  treated  as  an  eigenvalue  or  as  a 
mechul  as  described  in  Ref.  [7].  The  latter  procedure  is  much  more  efficient 
than  the  former  procedure  and  is  used  here. 


To  solve  the  BF  model  equations,  we  first  introduce  the  stream  function 
it>(x,y)  as  in  Ref.  [7]  in  order  to  satisfy  the  continuity  equation.  With 

=  w  treated  as  mechul,  the  resulting  system  can  be  written  as  a  system  of 
four  first-order  equations: 


w'  =  0 
<Ji'  =  u 


+  u  3U  _u.  !i.  p  +  T. 
at  ax  ax  3 


3/2 


lk+  u  !£  =  2al  TU'  “T~“Tmax(GT)' 


(29a) 

(29b) 

(29c) 

(29d) 


We  again  center  Eqs.  (29a, b)  about  the  midpoint  (x. ,  t  ,  .  )  and  Eqs. 

i  n 

(29c, d)  about  the  midpoint  (x_.  ,  ,  t  ,  ,  tv-  .)  of  the  cube  shown  in  Fig.  1. 
This  procedure  yields  the  following  nonlinear  algebraic  equations: 


(30a) 

(30b) 


i,n  1  ,n 

Tj  j-1  . 

FT  i 


u^-u1’" 


u]_1  ,n  —  ui ,n 


1 — J— y  ”  -  («*>]$ 

,  (30c) 
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'k  Ut\  ..1,n  ,_i.n  .  .1-1, ns  .  ..1-l,n_1,n  *  j  .1-1  ,  i.n  .1-1. 

Vj-Js  +  °1  [ ua-«s  {Ti-h  Tj-*  >+Vs  Tj-*sJ  al|  - Rp  (*j-*  “  *j-»s 

J  /J-'l.n  _T1“1.n\l  u1.n_  1,n  L  r,  3/2^. n  .  3/2»1-l ,n 

♦  ♦  -t -".v .  a=i  4  2  firJH :  «  ia 

j-*s\  n*  /  ,  l,n  “  ,1-l.n 

\  J  /J  J  Lj-Js  Lj->s 


i.n  1,n 

j.  /r>  i  \  1  »n_1  ,n  j_  »h  J  J  ■  1  _  _ 

♦  (G  IjVj-l,  +  Gj-H  J-  hj  "4 


(30d) 


where  now 


a  a  - 

1  X1  ~  X1-l 


1  *n  Vl 


*  _  °1 
ai 


•k  .  in_ 

6n  55^ 


"3  *  -*(-3)^  -  <*'  +  2*„<“U’"  -  “U’"-1  - 

"4  *  V*3  -  V  -  26nA1  +  -  2\  -  <G'"“  -  <Gt')0?-^ 

.  .  ,1-l.n  _  1-1 

*1  J-%  3->5  j-V 

*2  ■  '  -  "jii’"'1)  - 

(^/2)l-r’  *  (t3'2)’:’-"-1 


"4  , 1  ,n-l  , l-i  ,n-l 

j-H  Lj-*s 

Again  the  system  given  by  Eqs.  (30)  Is  nonlinear  and  is  linearized  by 
using  Newton's  method.  This  procedure  gives  rise  to  the  following  form 
(2  <  j  <  J) 


<Wj  "  SWj-l  =  (r3>3 
h. 


5*j  ~6Vi  ~r  (6uj +  6uj-i)  =  (r4>j-i 


(31a) 

(31b) 


12 


<r3>j  '  0 


(r,)3  *  "3  -  I2enuj-h,  *  _  “1 

+  (u'>5:i'VH>-(T'W 
-  '£•">  *  “U’Vs' 

-V<<'  w*^  - 

#  -  <^a- 


(r2)j  -  "4  - 


+  6jVj-H  +  <WT')j->i 


For  j  ■  1,  we  use  the  boundary  conditions  given  by  Eqs.  (13),  (14)  and 
(15)  and  first  write  them  as: 

u,  «  *^2.5  In  +  5.2) 

3W, 


3*1  u1y1  3W1 
ax  w’i  T5T 

u,  3w 


9  U1  ""I  3  2 

T1  *  W1  +  yl  *^W  *  “*yl  ax  W1  "  p3yl 


After  we  write  the  difference  equations  and  linearize  the  resulting  nonlinear 
expressions  we  get 


5u,-j2.5(ln^.^> 


5w1  =  (r1)1 


(32a) 


y1(w1  -  E2)6u-,  +  (w1  +  w^34)6\(i1  +  [(<1^  -  Ei)  -  y1(u1  -  u^34l)«w1  s  (f^)-, 

(32b) 

6t^  +  g7«w^  *  (r3)^  (32c) 


where 


il 
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e,  •n-1 

E2  =  (wj)i_1>"  -  (wf)'-"-1  ♦  (wf)1-1*"-1 

C  =  -W1'n"1  +  Wi_1»n  -  W1'1’"*1 

+  fc4  wl  W1  W1 

y-j  2 

E5  =  f  [2. 5 In ( 1 . 0  +  yj)  +  5.1  -  (3.39yj  +  5.1)  exp  (-0.37y|)J  dyj 


g?  =  -2w1  [1  +  J  yi°iE5234l  2  yl8n  (w]  ^ 


u,  1234 


(r1)1  =  Wl[2.5  In  +  5.2]  -  U] 

/  \  _  1234,  r  x  1234,,  c  \ 

(r2^i  ~  (w*i  ^2)  w  Et  ) 


'i  ^ 


(r3),  -  (wf)234  -  ,f 4  -  4(P3);:Jy,  *  \  y,6„(^-),234E4  *  \  y^E^E, 
For  j  =  J,  we  use  the  usual  boundary  condition. 


UJ  *  ue 


which  in  its  linearized  form  is 


«Uj  *  (>"4)1  *  0  (33) 

The  equations  (31)  for  2  <_  j  <  J  and  the  boundary  conditions  given  by 
Eqs.  (32)  and  (33)  form  a  linear  system  which  is  solved  by  the  block- 
elimination  method  discussed  in  Ref.  [9]. 
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IV.  RESULTS  AND  DISCUSSION 


To  study  the  calculation  of  unsteady  turbulent  boundary  layers  with  and 

without  flow  reversal  we  have  considered  three  separate  test  cases.  The  first 

one  has  an  external  velocity  distribution  of  the  form 

ug  =  1  -  a(x  -  x2)(t2  -  t3)  0  <  x  <  1 ,  t>0  (34) 

where  a  Is  a  positive  constant.  The  same  velocity  distribution  was  recently 
used  by  Cebeci  [1]  for  laminar  flows  in  order  to  study  the  computation  of 
unsteady  laminar  flows  with  flow  reversal  using  the  solution  procedure  described 
in  the  previous  section  and  to  see  whether  there  is  a  singularity  associated 
with  such  flows. 

In  performing  calculations  for  this  case  and  for  the  others  considered 
here,  care  must  be  taken  in  generating  the  initial  conditions  in  the  (t,y) 
and  (x,y)  planes  at  some  distance,  say  x  =  xQ.  For  a  laminar  flow  if 
xQ  =  0,  the  initial  velocity  profile  for  the  velocity  distribution  given  by 
Eq.  (34)  can  be  taken  as  Blasius  and  there  is  no  difficulty  about  computing 
the  solution  in  x  >  0  since  the  initial  boundary  layer  is  of  zero  thickness. 

If  xQ  1  0,  we  can  take 

ue  =  1  -  a(xQ  -  x2)(t2  -  t3)  0  <  x  <  xQ 

but  then  at  x  =  xQ  there  is  a  discontinuity  in  the  pressure  gradient.  Since 
it  acts  on  an  already-established  boundary  layer,  the  initial  response  is 
inviscid  leading  formally  to  a  velocity  slip  and  hence  a  subboundary  layer  at 
the  wall.  The  treatment  of  the  boundary  layer  is  then  rather  subtle  (see  Ref. 
[10])  but  if  we  are  not  too  concerned  with  the  details  of  the  solution  near 
x  =  xQ,  which  is  the  case  here,  a  convenient  procedure  would  be  to  write 
Eq.  (34)  as 

ue  -  1  -  aF[(x  -  xQ)/a](x  -  x2)(t2  -  t3)  (35) 

where  F  is  a  smooth  function  which  vanishes  if  x  <  xQ  and  is  unity  if 

x  —  xQ  >  a.  For  example,  we  can  take  F(s)  =  sin(irs/2)  0  <  s  <  1,  and 

a  =  0.06  with  ten  stations  between  x  *  x.  and  x  =  x  +  a.  A  similar  dif- 

0  o 

flculty  would  occur  at  t  =  0  if  Xr  in  Eq.  (34)  were  replaced  by  t  since 
the  boundary  layer  is  well  established  at  t  =  0. 
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Figures  3  and  4  show  the  results  for  the  turbulent  flow  calculations  with 
the  CS  model  for  the  test  case  given  by  Eq.  (35)  with  a  =  40  and  a  unit 
Reynolds  number  uQ/v  *  2.2  x  106/m.  The  results  shown  in  Fig.  3  were  obtained 
by  using  different  expressions  for  A;  those  shown  by  circles  were  obtained 
with  Eq.  (7),  and  those  shown  by  solid  lines  with  Eq.  (7)  written  as 

A  =26(1)  4  (36) 

p  max 

As  can  be  seen,  both  expressions  give  nearly  the  same  results. 

The  results  in  Fig.  4,  as  in  laminar  flows,  exhibit  no  signs  of  singular¬ 
ity  for  all  calculated  values  of  t.  This  is  in  contrast  to  the  findings  of 
Patel  and  Nash  [11].  Again,  as  in  laminar  flows  [3],  we  see  the  familiar 
rapid  thickening  of  the  boundary  layer  in  the  reversed  flow  region.  If  it 
had  not  been  for  this,  the  calculations  would  have  been  computed  for  greater 
values  of  t  than  those  considered  here. 


The  two  other  test  cases  considered  here  correspond  to  Cases  4  and  5, 
as  reported  by  Carr  [12].  Case  4  is  for  unsteady  Howarth  flow.  It  starts 
from  a  well-established  steady  flat-plate  flow,  on  which  a  linear  decelera¬ 
tion  of  ue  is  imposed  at  t  =  0.  The  external  velocity  distribution  is 
given  by 

ue  -  1  -o(x  —  1.24)t  1.24  <x  <4.69  (37) 

where  a  is  a  constant  equal  to  2.4/3.45  secern”"*.  The  flow  was  assumed  to 
be  steady  up  to  x  =  1.24m;  the  velocity  distribution  Eq.  (37)  was  then 
imposed  as  a  function  of  x  and  t.  This  test  case  differs  from  the  previous 
one  in  that,  once  the  flow  separates,  it  does  not  reattach.  For  this  reason, 
the  calculations  can  only  be  continued  as  far  as  the  station  where  the  flow 
reversal  first  occurs.  The  Initial  velocity  profile  at  x  =  1.24  and  for 
all  time  correspond  to  a  flat-plate  profile  with  a  momentum  thickness  Reynolds 
number  (RQ)  of  4860,  and  local  skin- friction  coefficient  c^  of  2.8  x  10  . 

As  In  the  previous  test  case,  we  Introduce  a  function  F^  so  that  at 
x  =  1.24,  dug/dx  =  0.  Since  we  also  want  the  solutions  at  t  =  0  to  cor¬ 
respond  to  steady-state  solutions,  we  Introduce  another  function  F2  In  order 
to  set  3ufi/3t  *  0.  With  these  functions,  Eq.  (37)  then  becomes 


Figure  3 


Local  skin- friction  variation  with  x  for  various  values  of  z. 
Solid  lines  denote  the  calculations  made  by  (29)  and  circles  by  (8) 


Figure  4.  Variation  of  displacement  thickness  with  x  for  various  values 


where 


ug  =  1  -  aF^x  -  1.24)t 


(38) 


t-  _  _J_  IT  /X  1.24\  C  -  rJn  Tt  t 

Fi  • s1n  r  (~tn — }*  F2‘siniT^s 

Figures  5,  6  and  7  show  the  calculated  local  skin-friction  coefficient 
cf,  the  shape  factor  H  and  the  momentum  thickness  Reynolds  number  R0  for 
this  test  case.  The  calculations  were  done  by  using  both  CS  and  BF  models; 
the  results  shown  by  solid  lines  refer  to  the  predictions  of  the  CS  model  and 
those  shown  by  circles  refer  to  the  predictions  of  the  BF  model. 

As  seen  from  these  three  figures,  there  is  essentially  no  difference 
between  the  predictions  of  both  models.  Although  there  is  some  discrepancy 
in  the  shape  factor  predictions,  this  does  not  seem  to  be  too  significant. 


t 


x 

Figure  5.  Computed  local  skin-friction  distribution  for  test  case  4. 
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According  to  the  predictions  of  the  CS  model,  which  also  has  the  capabil¬ 
ity  of  predicting  unsteady  boundary  layers  with  flow  reversal,  the  wall  shear 
vanishes  first  around  t  0.22,  x  =  4.69.  Since  the  computation  of  boundary 
layers  for  values  of  x  in  the  range  1 .24  _<  x  _<  4.6^  for  t  >  0.22  depends 
on  the  specification  of  a  velocity  profile  at  x  =  4.69,  we  generate  such  a 
profile  by  assuming  it  is  given  by  the  extrapolation  of  two  velocity  profiles 
computed  for  x  <  4.69.  This  procedure  in  which  the  extrapolated  station  serves 
as  a  downstream  boundary  condition,  allows  the  calculations  to  be  continued  in 
the  negative  wall  shear  region  as  shown  in  Fig.  8. 

The  third  case  considered  in  our  study  corresponds  to  Case  5  in  ref.  [12], 
which  in  a  way  resembles  the  external  velocity  distribution  in  Eq.  (34).  It 
is  given  by 

ue  =  1  +  )A2  +  (Bt)2  [£-S0]2}2-  [A2  +  (BCQt)2]^  (39) 

where  A  =  0.05,  B  =  3.4  sec  ,  ?  =  (x  —  1.24)/3.45  and  the  range  of  x 
values  are  limited  to  1.24  <_  x  £  4.69.  As  before,  the  initial  velocity  pro¬ 
files  at  x  =  1.24  for  all  t  correspond  to  a  steady  flat-plate  flow  with 
Re  =  4860,  Cf  =  2.89  x  10~3.  We  again  modify  Eq.  (39)  to  avoid  the  dis¬ 
continuity  in  the  pressure  gradient.  This  time  we  multiply  the  right-hand 
side  of  Eq.  (39)  by  F-|  used  in  Eq.  (38). 


Figure  8.  Variation  of  wall-shear  parameter  f^  with  distance  as  a  function 
of  time  for  test  case  4. 
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Figures  9  and  10  show  the  calculated  local  skin-friction  coefficient  cf 
and  the  momentum  thickness  Reynolds  number  RQ  for  this  test  case.  Again 
we  present  the  predictions  of  both  turbulence  models.  Figure  11  shows  the 
calculated  velocity  profiles  for  several  t  and  x-stations.  As  is  seen 
from  these  figures,  the  predictions  of  both  turbulence  models  are  the  same 
for  all  practical  purposes. 

Figure  12  shows  the  variation  of  wall  shear  parameter  f"  as  a  function 
of  x  and  t,  and  Figure  13  shows  the  calculated  velocity  profiles,  includ¬ 
ing  the  regions  in  which  there  is  flow  reversal  across  the  boundary  layer. 

These  computations  which  are  done  by  using  the  CS  model  provide  confirmation 
of  the  general  trends  in  test  case  4,  namely  that  as  in  laminar  flows, 
the  unsteady  turbulent  boundary  layers  thicken  rapidly  with  increasing  flow 
reversal.  A  new  feature  however  is  the  dip  in  the  graphs  of  f"  near 
x  =  2.5  which  develops  as  t  increases  towards  0.40.  It  is  possible  that  a 
singularity  occurs  in  the  solution  at  a  later  time  as  many  authors  have  sug¬ 
gested  is  the  case  for  laminar  boundary  layers.  The  most  cogent  argument  in 
favor  of  this  phenomenon  has  been  advanced  by  Shen  [13J  but  we  note  that 
the  most  definite  sign  of  its  occurrence  appeared  in  his  graphs  of  displacement 
thickness  which  showed  spikey  characteristics.  Here  the  displacement  thickness 
seems  to  be  fairly  smooth  but  the  skin  friction  becomes  spikey. 
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V.  CONCLUDING  REMARKS 


Based  on  the  studies  conducted  in  this  report,  we  observe  that: 

1.  The  numerical  solution  of  unsteady  laminar  and  turbulent  boundary 
layers  Including  the  flow  reversal  across  the  layer  can  be  obtained 
quite  satisfactorily  for  a  given  pressure  distribution.  A  combina¬ 
tion  of  both  regular  and  zig-zag  box  schemes  are  shown  to  yield 
accurate  results  for  unsteady  boundary  layers. 

2.  Whether  the  unsteady  boundary-layer  equations  for  laminar  and  turbu¬ 
lent  flows  are  singular  for  a  given  pressure  distribution  still 
remains  to  be  investigated.  The  results  for  test  case  5  indicate 
that  at  large  times  there  is  a  puzzling  "kink"  in  the  wall  shear 
parameter,  f^;  this  may  be  due  to  a  singularity  or  it  may  be  due 

to  a  numerical  problem.  Recent  studies  conducted  by  Cebeci  [14]  and 
van  Dommelen  and  Shen  [15]  for  a  circular  cylinder  started 
impulsively  from  rest  indicate  that  at  large  times,  t  =  1.25  or 
more,  there  appears  to  be  a  singularity  in  s*  around  $  =  120°. 
However,  these  calculations  do  not  indicate  any  puzzling  behavior 
in  thewall  shear  parameter  near  "singularity;"  the  f"-values 

W 

are  smooth  and  well  behaved  for  these  and  larger  times.  On  the 
other  hand,  examining  the  6*-results  for  test  case  5,  we  find  that 
while  there  is  an  abnormal  behavior  in  f"  at  large  times,  the 
corresponding  fi*-values  are  smooth  and  well  behaved,  a  trend  which 
is  opposite  to  that  for  a  circular  cylinder. 

3.  A  comparison  of  the  predictions  of  two  turbulence  models,  namely, 

CS  and  BF  models  Indicate  that  for  attached  flows,  both  models 
yield  almost  Identical  results.  This  is  also  true  for  flows 
which  are  sufficiently  strong  in  pressure  gradient  to  cause  flow 
reversal  across  the  layer. 
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VII.  DESCRIPTION  OF  THE  COMPUTER  PROGRAM  WHICH  USES  THE  CS  MODEL 

Input 

Essentially  the  Input  to  the  computer  program  consists  of  four  types  of 
cards.  Card  1  contains  the  title  of  the  flow  problem  under  consideration. 

Card  2  requires  the  following  Information  to  be  specified. 

i 

i 

NXT  Total  number  of  t-stations  to  be  calculated 

NZT  Total  number  of  x-stations  to  be  calculated 

NTR  x-station  where  transition  begins.  If  the  Initial  velocity 

profile  is  for  turbulent  flows,  then  NTR=1.  If  flow  is  all 
laminar,  set  NTR>NZT. 

IBDY  Specifies  whether  the  flow  at  x=0  starts  as  a  flat-plate 

flow  or  as  a  stagnation-point  flow. 

=1  flat-plate  flow 
=2  stagnation-point  flow 
RL  Free-stream  Reynolds  number,  u  L/v. 

oo 

IPRNT  Controls  the  print  output 

=1  prints  out  only  the  boundary- layer  parameters  6*,e,  cf, 

R5*.  R0,  H  and  external  velocity  distribution. 

=2  prints  out  profiles  as  well  as  the  boundary -layer 
parameters  and  external  velocity  field. 

DETA(l)  and  VGP  are  the  nonuniform  grid  parameters  that  control  the  spacing 
across  the  layer.  The  grid  used  in  this  report  is  a  reometric  progression  with 
the  property  that  the  ratio  of  lengths  of  any  two  adjacent  intervals  is  a 
constant;  that  is.  An-  =  KAn,_ r  The  distance  to  the  j-th  line  is  given  by 
the  following  formula: 

An,  =  h1 (xj  -  1)/(K  -  1)  K  >  1 

J 

There  are  two  parameters  in  this  equation:  h^  ,  the  length  of  the  first  step, 
and  K,  the  ratio  of  two  successive  steps.  The  total  number  of  points  J 
can  be  calculated  from  the  following  formula: 

lntl  +  (K-lKng/^  )] 
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In  the  computer  program  which  embodies  the  present  solution  method,  h^  and 
K  are  chosen  with  typical  values,  for  moderate  Reynolds  numbers,  of  0.01  and 
1.3,  respectively.  In  general,  approximately  50  grid  nodes  across  the  boundary 
layer  are  sufficient  to  represent  laminar  and  turbulent  boundary-layer  flows. 
The  chosen  values  of  h1  and  K  must  be  such  that  the  formula  which  generates 
the  number  og  rid  nodes  according  to  a  given  or  estimated  nfi,  l.e.  Eq.  (  ) 
does  not  allow  J  to  exceed  101.  Figure  14  Is  provided,  therefore,  to  provide 
guidance  In  the  selection  of  J. 


Figure  14.  Variation  of  K  with  for  different  -  values. 


CF  and  RTH  are  the  local  c^  and  Rg  values  which  are  used  to  start  the 
turbulent  flow  calculations.  The  initial  velocity  profile  is  generated  t>> 
using  the  formulas  proposed  by  Granville  (see  ref.  9) 


a  —  (In  y+  +  c  +  n(l  -  coswn)  ♦  -  n^) 

U  K 
T 


(40) 


From  Eq.  (40)  and  from  the  definitions  of  6*  and  e,  it  can  be  shown  that 

t*  }  U  -  u  u  u  ,, 

0 


t  e 
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+  T55-  ”  77  “  °*122925n]  (42) 

From  Eq.  (42),  taking  Si (n)  *  1.8519,  we  can  also  write 

(1.9123016  +  3.05603  +  1.5JI2) 

Evaluating  Eq.  (40)  at  n  =  1 ,  we  get 


For  a  given  value  of  cf  and  R0,  we  can  solve  Eqs.  (42)  and  (43)  for  6 
and  n  and  then  substitute  them  into  Eqs.  (40)  and  (41),  thus  obtaining 
u-profi le. 

Cards  3  and  4  read  in  the  t  and  x  stations,  respectively.  The 
present  computer  program  specifies  the  external  velocity  distribution  by  a 
formula  and  computes  the  dimensionless  pressure  gradient  parameters  analytically 
as  is  shown  in  the  listing  which  follows  this  section.  The  test  case  in  the 
computer  listing  is  for  case  4  of  ref.  (12). 

Output 

Depending  on  the  IPRNT,  the  computer  program  prints  out  the  profiles 
f,  f' ,  f"  and  b  as  a  function  of  the  similarity  variable  n  and  grid 
parameter  j  together  with  a  parameter  KALC(J).  Here  KAIC(J)  =  0  when 
we  use  the  standard  box  and  =  1  when  we  use  the  zig-zag  box. 

ETA  n 

F  f' 

U  f ' 

V  f" 

B  b  (=1  +  e*)  equals  1.0  for  laminar  flows 

The  output  also  Includes  displacement  thickness  6*,  momentum  thickness  e, 
local  skin-friction  coefficient  c^,  Reynolds  numbers  based  on  6*  and  e, 
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that  is,  R6*,  Rq  and  shape  factor  H.  The  definition  of  these  parameters 
and  their  computer  notation  is 


OELSTR 


6*  =  ^  (1  -  u/ug)dy 


THETA 

CF 

RDELST 

RTHETA 

RZ 

H 


o  =  f  u/uJl  -  u/u  )dy 
0  e  e 

cf  =  2tw/pUo 


R*  =  «*u  /v 
6  0 


Re  -  eu0/v 

Rz  =  2U0/v 


6*/e 


In  terms  of  transformed  variables,  6*,  e  and  c^  can  be  written  as 

6*  =  —  [n  -  f  /f  ’  3 

»  »  =D 

2 


6  =  f  it  (1  —  rr)dr 

AT  '  fe  fe 
x  0 


f 


*  2 


fM 

w 

z 


12 


C'-'MMrig/BiCn  / 

i 

CIMMPN/plC 1/ 

1 

r dmmpn/plCP  / 
l 

ppmmpn/FONT/ 


NXT,NZT.NX,NZ.NP,NTP,I  TKA X ♦  I  BUY  ,  NP7  ,  IZIG,  HUPP.  ET  Ac, 
VGP.AdOl)  .E’AllOll  .rFTAClOl) 

X(  1011  ,Z(81)  ,U0<81>  ,PZ(B1)  ,P1  Ol)  ,P2  (8  1 1,«>3 (101,81*. 
UP  1 101 »  81 1 

OFLVI 101)  ,F  1101 ,81,2)  ♦UUOl  ,81 ,21  .V  ( 10 !  ,81,2). 

8< 101, 81,2) 

TPPNT 


pOMMPM/jntd/  N  .  I  F  X ( 81) 

common /OA VP/  lTlMP,NXTn,NZT0,X0»101)  *20(81)  •  CCS(81)*H*‘(81)*P7HFTA 

roMMOM/pir^/  pp <  ioi » ,ur  ( ion  ,vc  (ioi)  ,bc non 

OTmfmS.JDM  FF(  101,3)  *UU(101  .3)  ,VV(101  ,3)  ,BB(10l.3> 
n TYCNS  ION  0  THC  TA ( 81 ) 


e c»r( 5, 101)  NXPCNT 

M  =0 

P^AX  =  10 
itimc  =  o 
5  MX  =1 

M7  *  1 

M  =0 

T~TM=  a  T-IMe  +  l 

CALL  TMPUT 

JC(  T T  TM F  .0-‘T.  1  )  GP  TP  f 

CALL  IVP 

nr-  io 

t  OP  7  J  *  1.M0  7 
M  J.  1 .  2 )*cc ( J  *  ? ) 
n<  J ,  1,2)=UU(J,3) 

V(  ).  l.?)»VV(J,3) 

°  *  J  «  l.  2 ) sPB (  J  *  ? ) 

7  -pn)~!MUc 

10  IC  (  ICC  v  ,LT.  2  .«rn.  TTIMP  .cr.  2 )  WPITF(6,110)  NX ,  N7  .  X  (  MX  J  . 
1  7  (  V  7  ) 

! c  (  rilfo  ,c.o.  0  l  nn  ’n  12 
T c  »  MX  .  rQ  .  1  .AMO.  M  7  .EO.  1  1  TO  TO  QO 

1  5  m\ (T  Tf; y p 

I F  (  M7  .-0.  1  .AMP.  MX  .AT.  J  )  C-C  T0  OO 

T*  =0 

’0OrV  =  0 
20  TT  a  !TM 

!cl  >T  .1  c  .  X)  GO  -n  70 

IOC  )  M  7 

o.p  Trj  ©0 

■»0  IF(»”  .c-c  .  MC)  CALL  PTPy 
"  (  MX  .0* .  1  )  r,p  to  5C 

c  »(  l  Tf  PM  7 

c,p  ”0  to 

so  csu  fpfcr,|  T  T  ) 
to  c  > ll  SCI  V7«  TT) 

CALL  vpp  tm 

-  pm  FrK  cpq  C'-’M  Vpo  OPNjr  c 

TFM  7  .cc  .  MTP  )  00  70  70 
'•—MV  'MAC  C|  oy 

»c(  are(cclv(  n  )  ,c,t.  c.oocnc-c  -c  20 
or  -p  po 

r  __vi|c  oiil  cmt  CL°W 

70  vFf  ARCfCFLVt  1)  )  ,L  ?.ec-0A  .pc.  APS  (DFLV  (1 )  / 
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i  ivii/,7,?)*o.5*nFivam  .it.  0.02)  go  to  so 

GO  Tp  20 

flO  TFfwo  .cO.  NPT)  GO  -P  c0 

!C(«RSflUMB-?.N7,2)/U(NP,N7,2|-1.0)  .  LT.  0.0015 
l  .AND.  AP<< V)NP .NZ.21)  .Lc.  0.001)  GC  T0  90 
TPMG*»nw.0c  .2)  GpTP  *0 
TGPfW  *  TG0  0W+-1 
WO I-=( 6. 120) 

IL  =2 
CALL  GROWTH)  IL  ) 
cr>  -n  ?(i 

c0  "ALL  ril'rDUT 


T  r 

( 

NX  ,GP 

.  N  XPO  N  T 

)  !PRNT  =  0 

T  p 

( 

IT!MC 

.9’.  2 

) 

GO  10 

T  c 

f 

JTTMC 

.50.  2 

• 

ANP.  NZ  .  F  0.  3  )  GC  ’0  <=6 

t  c 

( 

I  T  I  M  P 

.  Nc  .  1 

) 

GO  TO  10 

I^(  M  .l”.  3)  GO  TO  92 
nn  ci  j=1,mpt 
CP( J , 3)=c< J . 3, 2) 

HIM  J  f  3)=U<  J.  3,  2) 

VV(  J  .  3)=V(  J  .  3,  2) 
oo(  J  ,  3)=P  (  J  ,  3.  2) 
o 1  0°N” TMUr 

Gr  T)  c 

c?  M  =  M  ♦  l 

0*3  c  j  =  1  ,  K|  p  r 
oo(J,Hl  =  =(J,M,?) 

NIIU.HI  =  U(J.M,2) 

WV(J.M)  =  V<  J.M. ?) 

PO(J.M)  =  P( J,M.2) 

Ol  CniyJT  T K'  1  I c 

JC  I  M  ,r.  3  I  r.p  rr  10 
NN7  =  2 

'‘ALL  cppsnc  <  MPT.  oo,  UL.  VV,  n  B  »  NN'Z  ) 

'I  7  =  7 

r,n  -p  io 
9f>  **  =1 

r*p  97  J  *  1 . N°  T 
cc< j .m )  =  cr ( j  ) 

HIMJ.M)  =  Ur(J) 

vv)  j . w )  =  vc«j) 

97  PP(J.M)  =  °PCJ) 
nr  q  (J  «  =  1,7 

M  a  M  +  l 

9  P  J  =  l.MPT 
crU,«)  =  -<J,K,2) 

DIM  J.M  )  =  U(  J.K  ,  2) 

VV( J . M  )  =  Vt J.K , 2) 

9B  PPf J.M)  =  P  f J.K , 2) 

V  7  =1 

rtu  orpcpc  (  MOT,  FC,  uu,  VV  .  p B  ,  MZ  ) 

TTT^C  =  TTIMC+1 

TP  (  TPP'I”  ,LT.  2  .AND.  T^IMf  .GT.  1)  WCTTC(6,U0)  MX ,  MZ  ,  X  (  NX  )  . 
!  Z  ( v  7  ) 

G  ALL  PHTPUT 

34 


i 


GO  to  10 


100  lHOf  16Xf3?HIT^P^T10NS  FXCFEDFO  !  T  MAX  AT  N  7  =  -»  T  ^  1 

101  coom  4T(  I  5) 

110  (1H0.AHNX  =,I3,5X,4HVZ  =  ,I?,*X,3HX  = ,F10. 5 ,5X ,3H7  =,*r10.SJ 

120  cOPM  AT(  lH0»2Xt  * Q  DUN  PA  c  Y  LA  VF  P  HAS  GPOWN'l 

121  cpQMar  (  1H0. 5X , 5HN  X  =  .  I  3 , 5  X  ,  8H  X  <  NX1  =  «P1 0. F  »5X . 5HX I  =  .MO. 5 f 
1  1H0*  CX  ,  ?H  J  .5X.5HZ  (J1.5X.8HV  (  UA  LI)  ,1  0  X  ,2HC  c  , 

?  13X,  1HH.  12X,  fcH^TH=TA  f  13X,2HL'E  ,8X,6HFXTRAP  /  ) 

122  PhdmaT  (  1H  ,  5X,  I  3,  2  X ,M  1 .  fc  ,  5  f  2 X , E 1 3 . 61  .5X.I1) 

CM  p 
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«U«cnuTIMP  COFPG(IT) 

'“OMMON/RLCn  /  NXT,NZT.NX,NZ,NP,NTPtI  T  MA  X ,  I  P  DY  ,  NPT  ,  IZIG,  ITUPR,  ET  AE, 
1  VGP.A(lOl)  ,FTJ  (101)  .DETA(lOl) 

Covmqn/rlC 1/  X( 101) , Z(81 ) ,UC(81)  ,PZ(81)  ,P1 (81)  ,P2 (81) , P3( 101,8 1) , 
1  t)c  (  101 , 81 ) 

COMMON /BLCC /  Sl( 101 )  , S2( 1C1 )  , S3 (101 ) ,$4  1101)  ,S5 ( 101 ) ,S6( 10 1 ) , 

1  oi(  101)  ,R2(1C1)  .03(101) 

COMMON /PL CP /  PELV( 101) ,P (101 ,81,2) ,U(10l ,81 ,2) ,V(101,81,2), 

1  8(101,81,2) 
mMMnN/’CZC,/  KALC(lOl) 

COMMON/SOVc/  nzig.nzigs 


U(NP,NZ,2)  =  UE(N'X,NZ)  /  L'O(NZ) 

DOB  =  0 . 5*  (  UO  (  N  Z  )  +  UO  ( N  Z-  1 ) ) 

IP  (  IT  .GT.  1  .AND.  NZTGE  .EG.  1  )  GC  ’0  6 
N  Z  I G  =  0 
K ALC(  l)=0 
00  4  J=2,NP 

IIP  =  0.5*(U(J,NZ,2)  +  U(J-1,NZ,2)  ) 

I P( U 3  .GE.  0.0)  GO  TO  4 
N7 IG  =  1 
K  AL  C  (  1 ) ■  1 
r,n  to  6 
4  '■ONTiNiir 

6  OELX  =  X(MX)-X(NX-1) 

PCL7  =  Z(iMZ  W(NZ-l) 

7  P  =  0.5*(Z(NZ)  +  Z(NZ—1)) 

C EL  *  ZP/OELZ 
RPI.M  =  7  B  /(  DEL  X )  /UPR 
PPL2  =  0 . 5*  C  E  L 
fcL4  =  0. 2 5* CEL 
P  1H  =  0 . 5*  P  1  (  N  Z  ) 

P2H  *  0 . 5* P  2 ( N 7  ) 

PIP  *  0.5*(P  l(N7.)+Pl(NZ-l) ) 

P?°  =  0.5*(P2(NZ)+P2(NZ-i)  ) 

P 18H  =  0 . 3*p IB 

P  38  =  0.25*(P3(NX,NZ)  +  P3(NX-l,NZ)  +  P3(NX,NZ-l)  +  P3(NX-l,N7-l)) 

OIJPDZ  =  0.25*(U(NP,NZ,2)**2-U(NP,NZ-1,2)**2 

1  +LMNP  ,NZ,1)**2-U(NP,NZ-1  ,  1 )  **2 )  /  DELZ 

OUnX  =  0.5*(U(NP,NZ,2)-U(NP,NZ,i)+ 

I  U(NP,MZ-1, 21-U(NP,NZ-1 ,1) ) /DCLX 

°38  =  7.8*(PUSD7+DUDX/ljOB)-»-P2B*<U<NP,N7,2)**2+U(NP,NZ,l)**2 

1  +U(NP,NZ-1,2)**2+U(NP,MZ-1 ,1)**2)*0.25 

!  F  (  N  7.  .PQ.  N  Z  T )  GO  TP  10 
M  7  P  1  =  N'Z  +  1 

cl  =  Z(NZ)*<  (Z(NZPl)-Z(NZ) )/(Z(NZPl)-Z(N7-l)  )/0FL7  ) 

E  1 H  =  0  •  5*  e  1 

E  2  =  Z  ( N  Z  )  /PEL  X /UO  (N  Z ) 

E3  =  Z(N7)*0FLZ/(Z(NZP1)-Z<NZ-1))/(Z<NZ)-Z(NZP1)) 

5  38?  =  0 . 5* (  P3(NX,M7)  +  P3(NX-1,NZ)  > 

PUOX  =  (U(N0,NZ,2)-U(NP,NZ,D)  /PELX 
OUEr Z  =  0.5*( (U(NP,NZ.2)**2-U(NP,NZ-1 ,2)**2) *E1 
l  +  HHNP  ,NZ,  l)’**2-tMNP  ,NZP1 ,1)**2)  *F3) 

0387  =  Z(MZ)*(OUSPZ/Z(NZ)+DUDX/UC(NZ)K- 

1  P2(NZ)*(U(NP,NZ  ,2)**2tU(NP,NZ,1)**2)*0.  5 

10  PONT  INUE 

On  50  J=2,NP 


CR  »  0.5*(e(J,MZ,2>+MJ-l,N7,2M 

^VR  *  0. •>*(£<  J,NZ.2)*V(  J,NZ,21«-F  (J-l  ,NZ,2»*V(J-1  ,  NZ  .  2  U 
n«  =  0.5MUIJ  ,NZ,?)+U(  J-l  ,NZ,2I) 

USB  =  0.5*(U( J,N7. ?}**2+U( J-l ,NZ,2)**2) 

U*>’  =  0.*>*(U(  J.NZ-1 ,2>+U( J-l  ,NZ-1  ,2)» 

UR  4  =  0.*5*(U(  J,MZ,1)  +  U(J-1  ,NZ,1)  » 

VR  =  0.5*(V( J,MZ.?1+V( J-l , N  Z, 2 1 > 

nf=n  V  =  (R(  J,NZ,2)*V(  J,NZ,2»-R  (J-l  ,NZ,2I  *V(  J-l  ,NZ,2n/9~TS(  J-l  1 

CP4  =  0.5*(F(J,M7.1)+*:  (J-1,NZ,1H 

UCP4  a  0.5*(U(J  N7,l)**2+U(J-l,NZ,l)**2) 

TC(M7  .cc.  *ZT)  GO  TP  20 
!c  (  N 7 ! G  .e0.  1  )  GC  TO  3C 

TO  RVJ2  =  F(J,N7.11*V(J,NZ,1H-C(J,N7-1  ,1)*V(J,N7-1  ,l)  + 

1  C(J,N7-1,2)*V(J»NZ-1  ,2 ) 

=  VJ1  =  P(J-l,NZ.n*V(J-l,NZ.l)  +  F  (J-l  ,NZ-1  ,!)*V(  J-1,NZ-1,  n  + 

1  P(J-1,N7-1,21*V(J-1,NZ-1,2) 

=  c  !1  =  0.2*5*(  F  (  J,N7-1  ,1  )+c  (  J-l  ,NZ-1  ,1)  +  F  (  J,NZ-1 ,2)+F(  J-l.NZ-1,2)) 

=VP234=  0  .  *>*(  c  VJ  ?+c  VJ 1 1 

u  S  Q  2  =  0.5*(U( J,NZ-l,2»**2+l( J-l ,NZ-1 .21**2) 

m?p?  =  0.5*(U(  J,MZ-1,1)**2+U(  J-l  ,NZ *1  ,1)  **2) 

1J«- °  T  1  =  0 . 5*  (  USB  2+USB  3  I 

'JDK  1  =  0.2C*  (  u(  J,MZ-l,n  +  l)(  J-l  .N7-1  ,11+U(  J«NZ,1)  +U<  J-l,  77 , 1  n 

MDJ2  =  IJ(  J  ,NZ,  1)**2+U(  J  ,NZ-l  .1)  **2+U(  J,NZ-1  ,2)**2 
U'Jl  =  U(  J-1.NZ,1)**?+U(  J-l  ,WZ-1  ,1>**2+U(  J-l  ,NZ-1,2)**2 
ncR2’4=  O.B*(USJ2+lJSJn 

VJ1  =  V(  J-l.N'7  ,  1UV<  J-  1  ,NZ-l  ,1)+V(  J-l  ,NZ  — 1  ,2) 

V  J  2  =  V  (  J  •  7  ,  1>+V( J.NZ-1, 1)+V(  J.K'7-1  ,?) 

V  p  2  d  4  =  0.5*rVJ2-*-VJl) 

RVJ  1  *  c  (J-  1  »N  Z  •  1 ) *  V( J-  1  tN'T  ,1)+B(J-1  ,  NZ-l  ,11  #V ( J— 1  ,  N7  -  l,  1  1  + 

1  R<  J- l.NZ-1  ,  21*  V(  J-l  ,NZ-1  *2) 

BV.J2  =  P(J.MZ,1)*V(  J,N'7,ll+B(J, NZ-l, 11*V(J, NZ-l  ,n  + 

1  R (J,NZ-l,?l*V(J,NZ-l ,2) 

rvi  =  UcB4-T.C*USBT  1 

rvo  =  iJR2-  ?•  0*(JB  K  1 

rtf  a  =  cr  a-  2 .  0*RR  I  l 

-yp  =  (  °VJ2-nVJH  /pdta  (  J-l ) 

r-MPj  =  C**  l*Ccl-0.  5*0^1*  VR234*C  ^4+2. 0*PC  LM*CM? -CM8-P1B*pVP234 
1.  +  °  2R*IJRP  234-4.  0*P3R 

rrP-cIC.  IENTc  cpc  THF  D^r-tLAP  R  r X. 

M(J)  =  B(  J,N7  ,  2) /DCTA  (  J-l)  +  PlPH*c  (  J,NZ  ,2)+fEL4*(  PB  +  CM4) 

S2(J  1  =-0(J-l,N7,21  /PCTi  (J-n  +  PlRH*P  (  J-l  ,NZ,21  *C  <=  14  *  (  F  8  +  C  M4 1 
=  P  1BH*V<  J«N'Z,21+rRlA*  (VB+VE234) 

P  4( J  (  =  P 1PM*V( J- 1,N7 ,2) +GFL4* ( VB+VB234) 

5 (  J  )  =-(  P2P  ♦CEL  )*IJ<  J  ,MZ  .2J-BFLM 
«6(J  )  =-(  D2B+rEL  )*U<  J-l  ,NZ.7)-Rf  LV 

P2(jl  =  rN3J-(nPPBV*PlF*PVB-(P2B+r':L)  *USB+rcLT'MVB*PB+VB234«FB 
1  +'*M4*V«  )-BFLM*2.0*LP) 

«&'.'•(  J  1=0 
nn  to  AO 

30  cc  T  =  0.5*(F(  J,N?-1  .2H-MJ-1  ,NZ-1  .2)1 

V  d  2  =  O.B*(  V(  J.NZ-1 ,2)+V(J-l  ,NZ-1  ,2)) 

cvPA  =  0.c*(  c(  J,N7 ,1)*V(  J.MZ  ,1H-F  (  J-l  ,NZ  ,1)  *V  ( J-l  ,NZ,  1)  ) 

HPfr  =  0.5*(U(J.N7Pl,ll+l'(J-l  ,N7P1  ,1)) 
c  P  6  =  0. 5*(c< J ,NZP1 ,1) *F ( J-l ,NZP1  ,1) ) 

"RAF  =  0.2B<‘(U(J.MZ,1H-U(J-1,NZ,1)+U(J,NZP\  ,  l  H-U  (  J-l ,  N7  ®1 , 1 )) 
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VA4  fe  =  0.25*(V<  J,MZ,l)  +  V(J-l.NZ,n<-V(J,NZPl  .1)+V<  J-l,N?Pl,m 
0  PBV4  =  fe(J,N?,l)*V(J,NZ,l)-B(J-l,NZ,n*V(J-l.NZ,in/?='rA(J-ll 
r\  =  p  1H*  (  VB  2*FB2-UP2**2) +F3’M  UJB4-UP6)  *IJP46-VB46*(  P04-PBfcn 
r  2  =  OP P V4+P  1  ( N Z )*c  VP  4-P2  (N  Z)  ♦US84+2. 0*P3  BZ 

T3  =  -C2  +  ?.0*C1-2.0*F2*UB4 
r  ropPF  jr  IFNT<;  PfC  THP  ZTG-ZAG  PCX. 

Sl(J)  =  B(J.MZ,2)/r'B‘pA  (J-1)  +  P1H*P  (J,NZ,2)+F1H*(PB-PB2) 

S?fJ)  =  -B(  J-1.MZ.2) /PrTA  ( J-1)  +  P1H*P  (  J-!  ,MZ,2>+P1H*(  PB-PB2) 
c?{J)  =  P  1H*V(J,NZ,2)+P1H*(VB+VP2I 
P4»J1  =  o  1H*V<  J-1,NZ,2)+E1H*<VP+VB2) 

SC(J1  =  -P2(\'Z)*U(J*H7«2)-B2-Fl*l)(J»NZ»2) 

SMJ1  =  -P2(t'7)*U(J-l»N'  Z*2)-E  2— F  1*U  ( J-l  *NZ  »2  ) 
c  2(  J  )  =  C3-<0  =  =BV+PI(N71*CVP-P2(NZ1*USP-2.0*F2*IJB-E1*<U3**2 
1  -VP*Fe-VB2*FB+PB2*  VBU 

K  ALT  (  J  |  =  1 

AC  oi(J)  =  P(J-l,N7,2)-pU,NZ,2)+OPTM  J-1)*UP 
P  31  J-!  I -III  J-l.NZ  t2)-U<  J.NZ,2H-0F7ft  f 
BO  rCMT  IMUc 

TP  (  I"  .-0.  1  )  N7TGB  =  NZIG 
P3(K’P)*  0.0 
A  1(  1 )  *  0.0 
r?(l)  =  0.0 

c  cryp\i 

CV  p 
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8U8P0UTI\'C  ?DnY 

COMMON /BLCD/  NXT  «N7T .NX.NZ.NP »NTP , I TPA  X«  T BDY  »NPT . I  2  I G  * [tUkS, FTAP, 
1  VGP.AC101)  «PTAC10H  ,DPTA(lon 

rr>MMOM/RLC  1/  XC 101) ,ZC81) ,UCC81)  »PZIP1> ,P1 C81)  ,P2  C81) , P?C 101»P1), 
1  Uc ( 101  *  81 ) 

rnMMHN/BLCP/  DFL V( 1 01 ) ,F ( 101 , 81 ,2 ) ,UC 101 ,81 .2 ) .V ( 10 1 ,01 , 2 ) , 

T  PC  101 *81*2) 

roMvrN  /RLCO/  PL.CP.PTH 
n  jw  fm  $ inM  PPVCIOI) 
dimfmcipm  TB  (  101) 


GfiMTR  =  i.o 

TP  C  ITUOP  .NP.  0  )  GO  T  12 

Ur* T  *  0.0 

Ifl  *  1.0/UC(NX,NTP-1) 

OC  10  !  =N  TR  ,  N  Z 

IJO?  *  1.0/UPCNX.T) 

1JPI  =  UP !  +  ( UP  1+UP2  )*0.5*CZCIl-ZCI-ll) 

10  UP  1  =  Un  2 

r-r- i  =  R.35C-04*UF  (NX»MZ)**3  /  (IF  C  KX,NTR-l )  *Z  (  NT  o-l)  >  **1 .34 

cvpt  -  r,p*cl**o.6f*CZCMZ)-ZCNTP-l))*UPI 

TF  C  PXPT  p  .  10.0  )  GSM7P  «  1.0  -  FXPC-FXPT) 

12  CU'U  =  0.0 

FI  =  UC l.MZ, 2) /UC^P ,N Z,2)*fl. 0-UC1 ,NZ,2) /UCNP,NZ,2)  I 
m  13  J=2,mp 

c2  =  UCJ,N7,2)/UCMP,NZ,2)*(1.0-U(J«M?,2)/UCNO,NZf2)) 

*11*1  =  p|)M  l  +  (  F  l  +  F  2)*A  C  J) 

IP  FI  =  c2 

TH^i  =  SOP^CZCM?) /C»L*UOCK,Z))  >*St’Pl 
*  ’•  =  UC(NX.NZ  )  *  TWP  ’’A*  P  L 

TP  r  .IP.  42*5.0  )  r,r  TC  14 
TF  (  ”  .r,T.  6000.0  )  r,0  -*0  15 
X  P  T  *  O-/425.0-1.0 

dt  =  o.55*<  1. 0-c  XP  1-0  ?43*F0Pf  CXPI  )-2.PP*XP!)  ) 

PC  ’P  20 

14  PI  *  0.0 
Gn  Tn  20 

15  ° I  a  0.55 

20  TFLG  =  0 

0 ?  3  a  fqot  CHOCNZ)*ZCNZ)*dL) 

074  =  SO®  T  (  0  Z  2 ) 

VWflx  a  VC 1.NZ.2) 

OP  3C  J  =  2 . N° 

IFC  4RSC  VC  ).GT.CPSC  V»»A  XI )  VPAX=  VCJ,NZ,2) 

30  oON-INUr 

FOVr  *  0.0168*  C l. 55/C 1 . C+PT ) )*p Z2* (U (NP , NZ ,2 ) *FT A ( NP)-p(NP. NZ , 2) ) 
1*GA*tf 
J  =  1 

80  I  cf  T  CL  G  .=0.  1)  GO  TO  80 

PPLU8  =  (  P?(MZ)  /P  ZM*  (Uc  (NX,NZ)  /  UOCNZ)  )  **2*C1 .0/ A8S (VMAX )**1 .5) 

YPS  =  c  Z4*FTA (J)*  SOP  TCAPSCVMAX)*C1.0-11.8*POLUS))/26.0 
CL  =1.0 

IFCYPA  .IT.  4.0)  FI  =  q.C-FXPC-YCA)l**2 

FOVI  =  0.1fc*<*72*AB$CVIJ,NZ.?))*Fl*GAMTR*ETACJ)**2 

TCC FnV  I  .LT.  FnvO)  GO  TO  ICO 
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0 


7  eL  G  «  1 
<90  cov(  j  )*  CP>vn 
go  -n  no 

100  enV<  J  >*  5rV! 

IF»J  .1_ c  .  2 1  GC*  r0  no 

IFI  C0V(  J  )  .GT.cnvt  J-  1 )  J  GC^O  110 

cr>v(  j  )*  prv( j-i}*<fpv< j-n-pf'vu-r?)  »*vg» 

!  =  (  FOV<  J  )  .LT.  =  r:V'->l  GOTO  110 
-nv( J 1=  ?nvn 
T  PL  G  =1 

110  MJ.NZ.?)*  1.0  +  EnV(  J) 

-I  *  J  + 1 

TC(J  .IF.  M°)  GO  TP  £  0 
MOM |  =  NP-  1 

TR< 1)  =  0( 1,N7 (2) 

n«->  1  EC  J  =  2.N° 

TB(J  I  =  («< J- i,NZ, ?)+P ( J«N7,2) ) *0. 5 
15P  CPM^TNUc 

TO(VC)  s  TR|NO«l) 

no  17C  J=  ?,mom1 
R(J.N7,2>  =  ( -rR  (  J  )*TR  (  J+l  )  )*  0.  E 
1?0  roMT  IM»l= 

R  (  N  ®  ,  N  Z  *  ? ) =R ( NDM 1 ♦ N  Z  *21 

pptupn 

CN0 
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SU^niTTM*:  GSir 

roMMON/PLCO/  N XT ,N ZT,NX ,  NZ , NP  , NTP ,1  TMA  X ,!  PrV  ,NP’  ,  ! Z  I G»  I^UCB.  FT  AE» 
1  VGP  »A  (101)  »F  7/.  f  1 0?  )  *DE  TA  1101) 


! c(  (  VGP- !  .0 )  .LF.  O.CCU  GO  TC  10 

MP  *  /iL0GUETAF/rFTA(l))*IVGP-1.0)+1.01/ALCG(VGP»  ♦  1.000  1 
1)  =  FT&C*< VGP-1.0) /(VGP**<NP-l)-1.0) 

RP  Tn  20 

lo  md  =  p-'tp /of ta m  ♦  i.ccoi 

20  IC(NP  .1 c  .  101)  GO  TO  ->0 
W»!‘rc(fc,  50  ) 

STOP 

30  CT A (  1  )=  0.0 
no  AC  J  =2. 101 
PF-A ( J l=vGP*0F  TA ( J- 1 ) 

A(.l)  =  0. 5*PFTA  ( J-  1 1 

A.0  CTJIJ);  ct;  <j-  n+PCTA  U- l ) 

C  STyPK| 


*0  AT(  IHO.  *  MP  EXCEcrEO  lot  —  PFCGPAM  T FPMT NAT CD» ) 
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r 


SlIAonij-M^p  C,anwTH<Lll 


r 


roMMON/PL  Cr>  / 

1 

fowMOM/etCP  / 

1 


MXT»N7.T,NX,NZ,NP  ,NTB,I  TPAX,IPOY  ,NPT,IZIG,  mjPB,  ETAP, 
VGP,A(101)  ,ETA  a  Oil  «DE  TA  (101) 

"ELV( 101) ,F (101 ,61,2) ,U(101 ,81 ,2)  ,VI 101 ,81,2) , 

« ( 101,81,2) 


MP'’  =  VP 

MO  1  =  MD  4-  1 

M  P A  X  =  101 

!«=  (  LL  ,cO.  1  )  r,n  Tn  <35 

'IP  s  ND  ♦  2 

TP  (  NP  ,GT  •  NPT  )  NP=NPT 
O’PMAX  =  »'0 

«  0n  IOC  J  *  MP1,NPMAX 

r( J,MZ. 2) =  U(Npn,NZ ,2)* (CTA (J)-FTA (NPO) )+F (NP0,NZ,2) 
HI J.MZ, 2)*  U( N°0 ,N  Z , 2) 

VIJ.NZ, 2)  =  0.0 
B( J,MZ, ?)=  P (MPO.N7 ,2) 

100  Pr’N|T  INHE 

CMO 
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CURP  0 UT IN  F  I COM  7 

rriMMON/BLCO/  MXT,NZT«NX,N7,MP.MTptiTWAx,IPDY  .NP^.IZIG,  mi R B.  F”  AE . 

1  VGP  .  4  (  1 01 )  ,  E  TA  (101)  ♦DFTA(lOl) 

'•OMMPN/RLC  1/  X< 101) ,7(81) ,UO(81)  ,PZ<81)  ,P1  (8  1)  ,  P?  (8  1 )  ,  P3  (  10  1 , 8  1 ) . 

!  UP (101,81) 

COMMON /PL CP/  PPL V< 101) ,e (101 ,81,2)  ,U( 101  ,*1 ,2) ,V( 101,81, ?) , 

1  8(101,81,?) 

COMMON)  /rlcc  /  SU  101)  ,  S2  ( 1  Cl )  ,S3(  1011  ,S4<  101)  ,85  (101 ) «S  6(  10 1 ) , 

1  81(101), 82(101) ,R3(101> 


8  CL  =  0.0 

I c ( N Z  .CT.  1)  BEL  =  0.5*(Z(NZ)  +  7 (NZ-1 ))/(Z(N7)-Z(NZ-ll) 
o  IP  =  P  1(N7  )+BFL 
p  2°  =  p  2(  N  7  )  +BEL 

on  3C  J  =  2,mp 

r  off  IM  IT  ION  oc  AVERAGED  QUANTITIES 

CP  *  0.5*(P( J,NZ,2)+P (J-l ,NZ,2) ) 

UP  =  0.5*(tl(  J,NZ,2)+U( J-l  ,NZ,2)> 

VP  =  0.5*<V( J,N7,2)+V( J-l ,NZ,2)) 

EVn  =  0.5*<P( J,N2,2)*V(J,NZ,2)+E(J-1 ,N?  ,2) *V ( J-l , NZ , 2 ) ) 
l)cR  =  0.5*(U(  J,N7,2)**2+U(  J-l ,NZ, 21**2) 

OFOPV  =  ( P(J,N Z, 2)*V(J, NZ, 2) -B (J-l  ,NZ,2)*V(J-1  ,NZ ,2 ) ) / OET A( J -1 ) 

I P ( N 7  .GT.  1)  GO  TO  10 

o P8  =  0.0 

CUB  =  0.0 

rvp  *  o.o 

CRB  =  -P2(NZ) 

GO  T0  20 

10  CP 8  =  0.5*(e(J,NZ-l,2)+P(J-l,NZ-l  ,2)) 

CU«  =  0.5*(U(J,N7— 1,2)+U(J— 1 ,NZ— 1 ,2)) 

rVR  =  0.5*(V(  J, NZ-1, 2)+V(  J-l  ,MZ-1  ,2)) 

C  FVR  =  0.5*(c( J,N7-1,2)*V( J,NZ-1  ,2)+F( J-l ,NZ-1 ,?  )  *V  (  J-l ,  N7-1,  2) ) 
CIJ8R  =  0.5*(U(J,NZ—  1,2)**2+IMJ— 1  ,  NZ-1, 2)  **2 ) 

r  PEP  RV=  (B(  J, NZ-1,  ?)*V(J, NZ-1  ,2)-P(  J-l  ,  MZ-1 ,2  )  *V  (  J-l ,  NZ-1 , 2  >  )/ 

1  DFTA ( J- 1 ) 

CL  R  =  CDFR8V+Pl(MZ-i)*DFVP-P2(NZ-l)*CUSB+P3(NX,NZ-l) 

CRB  =  -P3(NX,N7)+BEL*  (CFVB-CUSB)-CLB 

r 

r  CnrFFlc  JFh'Tc  0F  THE  0  T  FFEF  r  NCFp  MOMENTUM  FQUATICN 

20  MU)  =  R(J,N7,21/0=TA  (J-1)+(P1P*F(J,NZ,2)-BEL*CCB)*0.5 

S?(J)  =  -B(  J-1,NZ,2)  /OF^A  (  J-1)  +  (P1P*F  (J-l  ,NZ,2)-BEL*CFB)*0.5 
S  3(  J  )  =  0.5*(P1P*V(  J,NZ,2)  +  BEL*C  VB) 

<4(J)  =  0.5*(P1P*V(  J-l,N?f2)  +  BEL*CVP) 

R8(  J  )  =  -P2P*U(  J,NZ,2) 

S  6(  J  )  »  -P  2P*U(  J-  1  ,M7 ,2) 
r 

C  OFF  IN  IT  IONS  OF  RJ 

R 1( J  )  *  F(J-l,NZ,2)-F(J,NZ,2)+OFTA(J-l)*UB 
R3(  J-1)=U<  J-  l.NZ,2)-U(  J,NZ,2)+DETA(  J-1)*VB 
0  2(  J  )  =  CRB-(0ERBV«-P1P*F  vB-P2P*USe-BFL*(CpB*VB-CVB*F8H 
30  CONTTNUc 
R  l(  1 )  =  0.0 
R  2(  1)  *  0.0 
p  3(  NP  )*  0.0 
PETURN 
FND 
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CijBcn^T  imc  jmpuT 

'■'jVMPn  /8l  C°  /  NXT,N7T,NX,MZ.NP.N',P,ITYAX,IPrv  t>'PT  ,!7!G»!Tl,fMr'’Ar. 
1  VO.U101I  (1011  ,rcTt(ll?1, 

roMUPM/m  r  i  /  X(  ion  , 71P1 >  .Un<8l)  ,°Z  111)  .PI  <81 )  ,°2  IP  II  ,  P31  10  1. 1  II. 

t  ijc  I  101.811 

-ommOM/pl^P/  0^1 VI 101) <101 .81,2)  ,U< 101  ,81  .2) , VI 10 1.8 1.21. 

I  11101,81,2) 

rrMyoK/RL^n/  pl.pp.pth 

Cnyy-IM  /PBMT/  r  PPM  T 

rpMMGM /CA vc /  T"IME ,NXT0,N7TC,XC(1011  ,  7018M 
IMmcm*!?*’  T  T  tl  e  (  20) 


°  I  =  ?.i4i5<n 

MP"  =  101 

r«4x  =  io 

~tac  _  8>0 

!«o  *  ? 

IN-e  «  2 
PI  =  1.0 

8a  =  2.4/1  33.0*’. 4M 
’=1  TMMF  .c0.  1  )G0  TP  60 
P X~=NXTP 
*"7  r=M?’ro 

or  40  i  =  1 ,  \>  x  r 

40  V(  I  ) *yn(  I  ) 

On  42  T  -  1 , N  7  T 
42  M  n=Zr<  !  ) 
pn  T)  po 
60  !TIJPR  =  0 

ocar  1  6.  270) 

,£4n  (  5,  260)  NXT,Nr,NTo,!RrY,?l.!PSM,DrMl),VGP.rF,PTh 

Tc  (  m-tr  ,~Q  ,  I  |  HUP"  *  1 

crao  (  r,  ?«=  c  i  ixm,  r=i,axri 
c  c/.n  (  2P)0)  (Mil,  T  =1  ,NZT) 

MX  to =M  XT 
'i  7  Tn  =MZ  T 
OP  70  T  =  1 » Ni  X 7 
70  XM  T  )=Xt  !  ) 

nn  72  I=l,a'7T 
7?  70m*7(I) 

7 1  =  Z<  1) 
nn  74  1=1,? 

74  MM  »  7  l-O.Ol*  (  3-  I) 

NX  T  =  l 
N7~=3 

«0  l>n  loo  !=1.MXT 

XI  I)  =  X (  I  1  *  ?3. 0 
100  C^MTIN'ic 

IM  !T’^  .c0.  1)  on  TP  1 2 c 
WOI'PU.  ??0 )  TITL* 

WO!*31  6.  340)  ►'XT.NZT.MTO  ,IP0> 

W°!TM6.  ’421  cL,  PFTA(l),  VO®.  CM  P  TH 

120  nr  1 4 C  T  =  i.NZT 
MC  1,  I  1  =  1.0 
0  31  1. T ) =0 .0 

ODTPOuPp  G°  A  P)  I FN  T  PAF^wcTpc  P  pp  STcar)Y  PTfi-pp 


44 


°?1  T  )  *  0.0 

Pll  !)  *  0.5*1  l.0+P2ini 
140  <jr{  If  *  1 .0 

tf<  !T!«F  .EO.  If  Bi=TUBN 
SAMPif  tf<t  CA$F  4 
0^  150  K  X  l ,M  x  T 
nr  15C  I  *  1.N7T 

cuwr  *  i.o 
','c  *0.0 
PtJMT  *  1.0 

*  c.c 

▼ c C  7  <  I  >  .r,c.  1.34)  GO  Tn  143 

=(iwr  x  «!M(P!/?.OM(  zm-1.24)  /O.ll) 

nc  =  d) /2. 0/0.  l*COMP?/2. 0*1  IZt!  1-1.24)  70.1)  ) 

143  TP(XCK)  .or.  1.08)  GO  Tn  145 
ci)v-  =  «IN1P  I /2.0*  (  X(K)  71.98)) 
r'T  =  p  I /l. 98/2.  0*C0S1PI  /2.0*  (X1K)  /I  .98)  ) 

145  iJrfK.T)  *  1.0-«A*X1K)*<7  <1)-1.24)*FUNC*FUNT 

ntfcrx  =(  -  PA*X(K  )  *F  UNC  -BA*X1K)*  (Z  (1  )  -1.24)  *DP  )  *F))NJT 

0tJPrT  =1  -BA*  (  7  (  I  t-  1 .24l*FUNC)*FUN'r*0T*(-BA*(ZU)-l.24)*X(K)*FUMC) 

P3(k  «  I  )  =?  1  I  )/U0(  !  )**2*(UF  <K,M*f'l'EDX*OUEDT) 

150  r0NT  INI,C 
160  W°  TT F  16,  32’)  mxT 

WRITE  1  6,  326)  (  XII),I=1,NXT  ) 

WRITE  (6,  324)  N7T 

WOI^c  (  6.  326)  1  Z  ( I  )  ,1=1, NZT  ) 

R  F  T(JR N 


760  POCMA^l  ETB.SCIO.O  ) 

270  crovAT(20A4) 

290  =npMAT( rc  10.0) 

722  R°R M A T  (  7//1H0, 27HTAPLF  pc  I K PLT  X  FPCM  1  TO  ,  J3  /  I 
324  cnoM  at  (  1H0, 27HTA  RLE  OF  INPUT  Z  FROM  1  TO  ,  13  /) 

326  C0RM6T  1 1H  ,  7Xf  1 2P 1 C. 5  ) 

330  =CRM4T<  1H0,  20A  41 

340  cor^at(  ///  1H0, 17H**  C  A  Sc  DA  TA /1H0 ,3  X  ,6HNXT  =  ,I3»14X ,  6HN’Zt  =  ,!7, 

1  14X.6HMTR  =,13/  4X,6HIBDY  =,I3  , 14X , 6HI Z I G  = , 1 3,  14X . 

2  6HI TUR  B  = , I  3 ) 

342  chomaTI  1H  .3X,6HRL  = , E 14. 6 ,3 X  ,6HDE TA 1= ,F 14. 6 , 3X ,6HVGP  =,F14.6/ 
l  4X.6HCF  =,F14. 6,3X,6HRTFFT= ,F14.6  ) 

350  FORMAT! 1H0, 3X, 6HB9  = ,F 1 4. 6 , 3 X , 6HBA  = ,F14.6 ,3X,6HBL  =,F14.6, 

1  ?X,  6HZ0  ®,F  14.  6 /) 

END 
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ciipc  '"Ht"  !v  c  '  V° 


r-<MOn/pLfn/  NXT.Mr^x.f/.tr.K’S.i’MU.lPDY  APT.!/!r..,’"M‘.F’  tK 

1  vr,p,»non,ni(101)  .f'FHliOll 

r-MWP^/Pi  ' w  X<  101 )  . 7  »  «1  )  ,ir <  8 l >  .p 7 (81 )  .PI  (4 I  )  ,  P? (0 1) . P3  I  10  1. B  l  ). 

’  'IMIOI.PH 

rA«',n*i/»!  ro/  r'eLV|  101)  . f  (101  «P1«?)  »  U  1 1  0 1  .PI  .2  )  .VtlOl.51.2). 

’  «  t  1 0 1  .  P 1 .  2 » 

■-PMMOV  /Ot  ro  /  0  L  .  r  r  .C  THE  T  1 

rn  MMPPj  /  etvc  /  !T!MP.NXTr,N?Tr.XC<10n,?r(«l> 

r  f«cNc  jo,,  un(  ion  ,rp  <  irt  i  ,rr,  ( i oi  i  .she  ac  (ion 

"a-*  o».OK.r/3.i415S2fE,c.41,2.0/ 


TC  (  T’-ltCB  .FO.  1  )  or  TO  P 

L  ftM  I  N|  *  C  P  C  O  F  T(_  c 


f.'ll  Gf  ! o 

=  l 

rTAMOo*  0 . 2e*c  ’’A  <KP  1 
c- an  i i  .*5/ct4  (^oi 
r'°  ?  J  =  l . Np 

x  r  Ta  <  j  )  /F  TA  (  mp  ) 
crte  ?  =  FTi«*» 2 

c  (  J  «  N  7  «  21  =  cTtMPO*cT^p2t(2.O-0.c*PTAP?) 
1 1 »  2)  =  0.  c*r  Tl»*  (  3.  0-c  Ta  P  2  » 

VIJ.N7,?)*  ctS  U1  6*  (  1  .  C-c  fAP2) 

«< J.N7, 21=  1.0 
3  rc MTINUC 
5  FTl|ON 


’•uopi)i_E‘:T  ppnpiLC 


p 
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C7f*7)=  U-(MX.M  7  )*  7  (K'71*°  t 
tno/  =  coc  T(0  j <M71  ) 

r  p°  2  =  0.^rc 

c0C  P"1  2  =  ?oo  T  ( C  cn  2 ) 

!frlM  =  SOCp°2 
«P  xr  c  2=  SOc  7*  POC  c0  2 
FO  =r  =  SOC  fn?/°K 
vy  *0.1 
AN  =  sorrnp/DK 

44  =  AN*AN 

CC1  *  -1.9l’30lf*«A+ll./12.**N 

rr  7  =  »M- •3.05603*7 1 

rC  3  *  -  l.  5*4  A 

rr  4  =  1.0/aN-C-4L0C-(^0rcr‘’*P'rHCTi» 

ftfi  =  CC1*0.  5*rr  2*rr  At-C.2c*CC3*rC4**2 

A42  =  0.5*C0r  2+CC  3<TC4I 

4  7  3  =  0.2  5*00  3 

Yt.nr,  =  ALnG(  VY) 

cpr  =  YY-14A  1*AA  2*Yl.rG+AA3*YL0G**2) 
ne  =  1 .0-114  2*2.  C*AA3*nrO)  /VY 
HYY  =  pcr/np 

YY  =  YY-PYY 

t  Ff  a  n  pf  Pyy  )  ,LT.  0.00C01)  G°  fr  20 
! T  =  ITM 


10)  GO  TP  10 


TM  1-  .1’. 

->0  rrN~*N1JC 

=»?c  *  0.5*lrr  4*4LrG<VY)> 

o  oc,  r4,  c  th:ti  /vv 

c’-t  r  *  c  nf  l  TA  7S0R  7 

fr(]TiMr  ,r,T.  1)  Wfi  I  Tf  (6, 1?0)  Plf  ,PDFLTA,FTAC 

cc^pqkx  CO.O/FCXPF? 

CAM  Gc  10 

’C  J  *  1 »  N° 

7c<c”Afj)  ,GT.  ccG°Ck>  GC  Tr  40 

in  r  p»|  t  j»n  jc 

M0 - G 1  *  J 
Gn  "*  *  0 
40  vc  e0  l  *  J  -  1 
vc =0  ?  *  J 

r;p  'it^  *t  CP, CCD  DQnci^f  « 

^  7  7  =  =  "M  N°  c  G  1  )  /CTAIMP) 

5  ?  =  e  -S  (NO  c  G  1 1*  sc  *r*  2 

c)  «  1  .  /pk*  (  A  L°G  f  c  2 ) *C  *P!  F*  (1  .-CCS  I  P I  *  /  Z )  >  ♦  7  Z  *7  Z  *  (  1  .  -  Z  Z  )  ) 

rv  =  O.pc 

”  =0 

7-  -V  FCC  x  c  lHv-*’M>|°2»ry, 

~c  x  9  I-C  7/f  1  .♦  (  P  2*p  V  )**2  ) 

-'-V  =  c  CC  /r,c  /r  v 

rY  =  c  Y*  (  1  .-or  Y) 

"  t  *  or,  Pr  y  )  .1  '.C.OOl  .GC.  JT.rM?)  GOTO  70 
*  ’  ’♦  1 

<-r 

7  ->  !-•-•<  -  rvitsi'n 

\j~v  -  cnr  ci-  7*  cc  x1*  c  ? 

"  er  I  =  1  * *'9  e  g  i 

II  l.v?,1)  x  ^r'n,,({T(v(lTfTM*CTf|J)l/'‘Y| 

o  "I  7)  r  v”«  /  I  1 . 0*  1I!TN-»**C ’•A  1  JM  **2) 

ir|  >  »  cr  j  .rr.  “  o  )  rr  Tr  ICC 

77  X  '  ~t  (  *JC  r  G  1  )  /c  TA  (  VP) 

'  =  i<v  e  c-  J  tNi  7 . 7) /cc  crj«_  <  AlfG  )  <C  KG  c?  *  =  T  A  <  N®  cGl)  )  *°I  c* 

’  f  { .-P'nf  e  1*  ■>n)+7  7»2J*  U.-77)) 

V  t  i  o  t  c  *p  t  /c  <  Np  I 

o  c  j  =►  9  c  r.  7 .  v0 

”  =  c  *'  <  J  )  /c  TA  (  *  o  1 

'  ~  r  7  V)  r  *  '"'MOT*//) 

-  T  N  *  *  >  r  =  frjP  7  1  1  . -cp  SANG**  2) 

<f  I .  *'7.  7)  =  c-  rc-,n*  |  f  Lrr,  |  cd  xf  c  ?*c'rA  I  J)  )  ♦C^Pl  c*(  l  .  -CT'  A^GI  4ZZ  *Z7  * 

1  (1.0-77)) 

<9  0  \M).N7.->)  r  cr  cnyK*  ,  1  . /c  fft  (  j|  ♦  VJ*'T  NANG+ 77/FTA  f^P)  *  )  7. -3  .  *77  )) 

r  f.  I.  "  III.  4  T  c  r  p  Q  c  t  c 
100  r~c*  =  7  ,*r  y**  7*  cc  xr  c  7 

ef  ) .  "7  .  ?>  *  0.0 
OP  1  10  J  =  ?,  vc  Cf,  1 

in  -(  I  .07.  7)  =  fOr  "o?*(c  TM  J)  y*A  -(AMI  TKTv*cta  (  J)  )  -1  .0  /  crcv* 

’  M.->G«  l.OMUT*  TM*E  TA  <  J)  )**?)  ) 

’  C  (  (  0  f-r  1  .CO  .  NP)  GP  fr  I  3C 
*’41  x  cT,(wc:r,l| 

C  |  x  o  ]C*Cf4  (MP  )  /f  T 

77  x  PTM7cTJC(B) 

C-  =  Cf»crr.  1.N7.7I  -  CCCO^K*  (F ’Al  *  «A  IPG  )'C)(CC?  *E’ A1 ) -1  .0*C*P1  C)-CJ* 

r  TV  (3  !*  7  7)  ♦C  T4  ,  ♦  7  7*  7  7*(  1 .  n.  -7  7/4.  )  ) 
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DO  120  J=MRf=G?,NP 
cTfij  =  c^MJ) 

7  7  =  FT<  J/ET&(NP) 

120  c(  J.N7,  ?)  =FC**CcP2K*  (FTA  j*  (A  IPG(  Sp  XCF  2  *FT  A  J)  -1  ,0*  O  PT  c )  -  PJ  *c  IN  (  P  I 
1  *Z7H-ETAJ*ZZ*ZZ*<1.  /3.-ZZ/4.  )» 

no  no  1 4 C  J  -  l.NP 

140  UIJ.MZ.2l  *  IJ(  J.NZ.2)  /U ( NP , NZ  ,21 
V<NP,N7, 2)  =  0.00010 
CALL  FDOY 
0  |=TIIPN 


no  cOPMAT(  1H0,  4H0TP  =  ,F  14.6nx,7HROP  L^A*  ,  F 1 4.  fr  ,3  X  ,  EH  ET  A  F=  ,  PI  4 .6 , 3*  , 
1  3HrF  =  ,<:l4.6,BX,flH9THCTAl=,Fl4.6) 
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supROUT  IN  =  OUTPUT 


common /pl  m  / 

1 

CO mmpn/PLC  1/ 

1. 

-ommom /pl C°  / 

1 

rn^MfiN  /7r,?0  / 


NX-r,N7TfNX,NZ,NP.NTp,ITMAX,IBDY,NPT.  !  Z  I  G ,  ITU  F  9 ,  ET  AF  . 
VGP.A(lOl) ♦FTA(lOl) »nFTA (1011 

X< 1011 ,7(81) ,110(81) ,RZ<81) ,P1 (81 ) ,P2 (fll) ,P3f 101.P1). 
UP ( 101 , fll) 

OELV< 101) ,c (1C1 ,81,2) ,U(101 ,81 ,2) ,V  ( 10 1 .8 1 , 2  > , 
fl< 101, fll. 2) 

KALC ( 1 C 1 ) 


common /pL  CD  /  RL.CF,CTH 


CPMMPN/PPNt/  TPONT 
rnv»MPN/TMTp/  M,1PX(  81) 

common  /  cj  \jc  /  l  TTME  ,MXT0.N7TC,  XT!  1C1  )  ,Z0<81)  ,P  =  S(fll),H<;(«l),F'rH^i 
n  tvi-MP  jcn  o  Tpic-rA  r  PI  )  ,o  THT  (  81  )  ,K'PK  (  81 ) 


r-.vr.jc  jr>M  F  =  ( 101,2)  ,UU< 1C1 ,2)  ,  VV  (101  ,2)  ,PB  (101  ,2) 
r'Tveyt  Trsj  (  101)*  «  l  01  )  ,  VN < 1  Cl >  ,8 N ( 1 01 ) 


*"7  ?  =  *'  7  T-  ? 

rcyroo*  Q 

T  c  X  (  N  Z  )  =  C 

7  -  (  .Lt.  3)  DC  T  c 

’  F  (  JPfMT  ,nc.  0  )  GO  T0  p 
WO!T=(t.  ?20  ) 

MO'*L  =  VD-1 

J  =  1 

2?0  ,  J,=  ”A  (J1  ,F  <  J.NZ  .2)  ,U(  J,N'7.2)  ,V(  J.NZ.2)  ,P(  J.N7  .  2), 

!  K/L^fJ) 

UP«;  *  MD--3 

WC  TTP(  (  .  230  )<J,STA<J),C(J,NZ,2),U<J,N?,?),V<J,N7,?),0(J,N7,2). 

:  KALC ( J) , JsNPM? .VP) 

c  ^o»r  t»'|jc 

MPK( N7  )=N° 

IP  (  7  .30.  1  .AMP.  N  to  .  N'F .  I  )  nr  to  90 

Cl  =  C0°T(  7  (  m  Z  )  /  (  0  l*Ur  ( M  7)  )  ) 

PPL  c  TP  =  C  (MP  )-c  (NP  ,KZ  ,2)  /l’(NP,N’7  .2  )  ) 

r  F  =2.0*V(  l.vz,?  I  /  (  SORT  (PL  *Ur(NZ)*Z  fMZ)  >  *(UF(  NX,  NZ  )/'jr(NZ  )  >**’) 

C  F  F  (  M  7  )  =  ^  c 

ODcLST=  Uc(MX,NZ)*pfl$tP#PL 
sum i  =o.o 

C1  =  IK  l.MZ,  2)  /UCN'P  .N'Z.2)*<1. O-LU  ,NZ,2)  /IKNP.NZ.2)  ) 
on  1C  J  *  2,  N!° 

c2  =  U(  J  .MZ  .  2) /U(NP  ,N  Z  ,?)*  (1.  O-IM  J.N  Z  .2)  /U(Np,N'Z  .  2)  ) 

CI'MI  =  SUMl  +  (Fl4-c2)*A(J) 

10  =1  =  c  2 

’■  hc~  a  =  r  i*  sum  i 

pthcta(m7l  =  Ur (NX,N7)*thpTA*RL 
=  CTHC’i(M?l 
H  =  PEL  S ’’R  /THP  TA 

MS(fZ )  =  H 

tr  (  IF  y  to  p  ,GT.  0  .ANP.  NZ  .cC.  ^7*  )  GO  ’n  109 
°0  rALL  GPnvi-H(  1) 

! c  (  I T IM E  ,LT,  3 )  Gr  TP  ICC 
I  F  <  IDCMT  ,GT.  1  )  GO  Tr  10c 

VIC  TTC  (  6,  242  )  DFLS-s  .THETA  ,  of  .PDF L ST  , P THE~ A < NZ )  ,H, 
l  Uc  (  N  X ,  M  7  ) 

100  Ic  (  M2  .FQ.  MZT  )  go  to  ICQ 
!c(  ITI*‘F  .L " .  3)  GO  Tn  140 
IF(  M7  .L  T.  MZ  2)  Gn  r p  14C 
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M  *  M  ♦  1 

m  120  J  =  l.MPT 
=  cfJ,N>  =  ef  JfNZ,2)  7UF(NX,N71 
MIJ(J,N|  =  U(  J,NZ,2) /LF <NX,NZ) 

\/V(J,Nl  =  V(  J,MZ,21  /(JC(NX,N71 
120  8Pf  J,NJ  =  «<J,NZ,21 

TF  (  NZ  .LT.  IWZT-).))  0  0  T0  140 
nr  124  J  =  l, NOT 

to  (UtJ.^7.2)  .  L T.  0.  1  OQ  *r  1?0 
12  4  rp^iiMijc 
00  "0  140 

130  07  1  =  7(m7T-d-71NZT-2) 
nZ  ?  -  ? ( N  7T l-Z ( M  ZT- 2) 
nr  134  J  =  1,NPT 
nc  =  Fc  (  J  , 2>-c=< J,  1) 

"IJ  =  UU(  J  .  2>-UU( J,11 
nV  =  VV( J . 21-VV( J,  11 
no  =  PP( J , 21-oB ( J,  1) 

FN(J  )  ~  -  r  (  J  «  1  )  +0  Z  2*nc  /O  Z1 
ijMJ)  =  UU(  J,  ll+nz2*OU/071 
VMJ)  =  VV(  J,  1)+OZ2*OV/OZ1 
°  M  (  J  1  =  88  (J.  1)+OZ2*DP/OZ1 
134  orN^TNUF 

icy-cp  =  1 

!  CX  (  N  Z  +  l )  3 1 
143  M  7  =  N7  +  1 

!Pf  !-r:t>  .1 c .  0)  GO  T0  148 
nr  142  J  =  l.NPT 
k  Al  0  (  J  1  =  2 

C<J,N7,2)  =  P\'t  J)*UF  (MX.N7) 

U( J  «  N  Z  «  2)  =  UN  I  J)*UF  (NX.NZ) 

V  (  J  ,  M  Z  ,  2  1  =  VN  {  J  )  *  Uc  ( N  X  ,  N  Z 1 
14?  «f J.M7, 21  =  “N( J) 

TC  (  tPcnt  ,MF «  0)  GO  "’0  5 
WCJ-  =  (f,  220  1 

f.;  p  ^  1  =  k  p  -  1 

WF  rTCl  6.  230  )(J.FTMJ),c«J,NZ,2)fU<J,NZ,2)tV(J,N7,2),P(J»N7,?K 
1  KALr ( Jl ,J=1 ,NPM1 ,3) 

J  =  No 

WPTTC(f,  230  )  J  ,FTA  (  J)  ,F  ( J  ,NZ  ,21  »ll<  J  ,NZ  ,2  1  ,V  (  J,  NZ,  2)  ,  B(  J.  N7 . 2  1, 

1  KALC(J) 

on  ro  c 

148  TF  (NX  ,GT.  1)  GO  T0  160 
r  TNI’Til  GlIPFr  FOP  NF  X T  F"A'r10N 
nr  150  J  =  l»M0T 
c( J,NZ, 2 ) =  C(J,NZ-1,?1 
iff  J.NZ,  21=  lit  J.NZ- 1  ,2) 
v  (  J  ,  N  7  ,  21=  V(J,NZ-1,2! 

1^0  MJ.N'7.21*  P(  J.NZ-1  ,2) 

TF  (  MX  .CQ.  I  )  8 F  TU°  M 
r 
r 

C  nc-rcpyjMr  MP  cop  7JG7&G  H  =  VP 
n  ( N 7  .c0.  NZT)  00  Tn  j?0 
fp  (MX  ,c0.  n  or  Tn  170 
7F  (NP  .L T .  N°K(N7+1)»  NP=NPK(N7+1) 

0  CT(|QNj 


SO 


lfcO  MP  =  MPK  INZ) 

IF  (N7  .CQ.  U  OPTUPM 

If  IMP  .IT.  K'PK(MZ'H)  NP*NPMNZ-1) 

Tc  <  I7IO  .FQ.  0)  GO  TO  17C 

TC  imp  .LT.  MP  KINZ+l))  NP=NPK  (  NZ+1) 

170  U9  =  Uc ( N X » M Z )  /  Uc ( N  X  ,  N  Z-*  1 ) 
nn  1PC  J  =  l  ,NP  T 
c(  J » NZ  •  2)  =  F(J«NZ—  l»2)*l)P 
U(J*NZ*2)«  U<J,NZ-1 .2>*UR 
V( J , NZ  *  2) s  Vt J,NZ~1,2)*UP 
B<  J.N7.2)=  MJ,NZ-lt2) 

180  CONTVUF 
o  c-ijON 

190  I F I  N  X  .FO.  MXT  ,/!N  0.  I ’I  ME  .FO.  1)  RF'rlJPN 

tp  [  icpn-  .\|c.  2  1  Gn  T  194 
1  XI  =  XINXI/37.0 

(  6,  280  )  NX,  X  IN  X )  ,  XI 

no  193  K  =  1,  MZ 

W5!TC  (  6,  260  1  K,  ZIK),  VII  ,K,2)  ,CP9(K)  (K»  ,RTHPTt(K)  , 

1  IF  (MX.K  J  ,IF  X(K) 

iq?  '"OM’TMMIP 

104  ’C|mX  .c0.  MXT  .ANO.  I  "I  MF  .FQ.  3)  STTP 
MX  =  NX+1 

N  =0 

M  7  =  1 

SHIFT. 

9°  21C  K  *  1  , M  Z  T 
or  2 CO  J  =  l,MPT 
f< J.K, 1>=  c(  J,K ,2) 

UI J,K,  1)=  U(J,K,21 
VI  J.K ,  1)  =  VI J.K ,2) 

200  8IJ.K,  11=  P(J,K,2J 
210  CPN’IMUF 
GO  T0  160 


220  ciOMflTf  JHO,  2X,  1HJ,  4  X , 3HC  Tt  ,10X,1HC  ,13X,1HU»1?X,1  WV  ,13X,lHp,8X. 

1  4HKALC) 

230  focm4t(  ih  ,  I3,r10. 5 , 4C 14. fc,I 6) 

2u->  cnBva,T(  1H0,  7H0ELS7P  =,E14.fc,3X»7H7HPTt  =  ,F14.  6  ,3X  ,7HCF  **E14.6  f 
\  IH  ,  7HC9PL  9T=fPl4.6,3X  ,7HP’HF-^t=  .F14.6  t 

c  IH  ,  7HH  =  ,F  Ii.6,3X,7HL'F  =,F14.6  /  ) 

280  e<  '  (1H0,3X,CHNX  =  ,  I  3 . 5  X  ,  RH  X  ( MX)  =  ,F  10, 8  ,  5X  .  5HX  I  =  ,F10.5/ 

1  1H0,  3H  J  ,5X,5HZ  (Jl.BX.PHV  I  WALL)  ,10X  ,2HfF  , 

2  13X,  1HH,  1 2X  ♦  6HR  'rHF  TA  ,13X,2HLF  »PX,6HFXTPAP  /  ) 

260  f^pmaT  ( 1m  ,8X,I3,2X,Fll.f ,5(2X,E13.6) ,8X,!1) 

c»  n 
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cttqcntjTTM?  CMrnTH 

rpMMnw  /»»(_  rn  /  NXT,nzt,NX.NZ»NP.NtR»!  tMAX,I  BDY ,NOT , 1 Z IG,  1tURB.  Ft  AE 

\  vr,o,A(ion  ,pta(iou  ,deta(iod 

rnMwr>N  /pl  r  i  /  X  (  101)  ,  Z  (  81)  ,UG  <  8 1)  .P  Z  (8 1 )  ,  PI  (8  II  ,  P2  (8  1 ) ,  03  (10 1 . 8  1 1 

:  uFdoi.en 

rnvMPW/PL  rp/  O-LVdCl)  ( 101  .81 t2)  .U<101  ,81.21 ,V( 10 1,8 1,2), 

i  mioi.81,2) 

p jmcmc  jnw  es»101l,U«tlCn,VR(101I.PS(101) 


NJPMJ  s  MP-  1 
MOM  2  =  MO- 2 

J  M  A  X  =  \ 

VMty  =  V( 1.NZ.2) 
r  1C  J  =  2, Mo 

T c  (  V<J,M?,?I  .LT.  VMAX  )  T-C  T  10 
WAX  =  V(J,M7,21 
JMAX  =  J 
10  COM‘r  ?MUr 

"UJ1  =  U<  JM?X+i,MZ,2)-U(JNAX,NZ,2> 
n/J  1  =  VlJMf  X+1,N7.2)-V<  JYAXtN7,2) 

JS  =  JMAX+2 

?c  J=J*.M° 

JJ  =  J 

0DJ2  =  U( J,NZ,2)-U< J-1  ,N2.2) 

"VJ2  =  V<  J,MZ.O)-V(  j-1  .NZ,2) 

U J  pc  on  =  UJ2*UJ1 
Vjonn  =  VJ  2*  VJ  1 

!<=  (  UJP  =  rD  .LT.  0.0  .OP.  vjpprn  .LT.  0.0  )  OC  T  C  30 
HHJ!  =  HJJ2 
f'VJl  =  n  vj  2 
20  nMTIVHiF 

20  !c  (  JJ  .c0.  MP  )  pcTtJpN 
nn  4C  J  =  JJ.MP 

PMJ)  =  0.5*(P(  J-1.M7.2H-F  <  J.NZ.2H 
UC|J)  =  0.  5*»Uf  J-1  ,N7,2)+Ut J.KZ.211 
VMJ)  -  0.5*(V<J-1  ,NZ,2)+V<J,A7,2> > 

»MJ)  =  0.5*(P(J-l,MZ.2)+R(J,f'Z,2n 

40  r<"M-r  tnuf 

no  *=;c  j=  j  j  *  m  p  m  i 

c(  J.N7, 2)  *0.  5*(  FS(  J)tf'(  Jd>  ) 

»)(  J.NZ.  ?l=0.e!*<US<  JI*U'I  J+1H 
V<  J.NZ,  2)=0.5*<  VS<  Jl+V5(  J*ll» 

B(  J.NZ.  2)  *0. 5*  (B2U)  TP  *<  J+IU 
•50  CrMTINUe 

V»o  =  -V»NP-l,NZ,2H-IUfWP.MZ.2)-U^PMl  .M7.211/ Af  NP1 
I-(ABSIVNP)  .LT.  APMV<NP.NZ,2I)  I  V<NP,M7.2>  *  VNP 
OC^UON 
cm  n 


r 


r 


«t!RPOUTIMF  CFPRPF  (  MPT ,  FF,  UU,  VV,  BB ,  K  ) 


OnvMnN  /pl  0  1/ 

1 

C^mwom/plCP  / 

1 

rn^MnN/poKiT/ 

roMMOM/pLC?/ 

ntwrfjc 

nj^FNCirM 
OTMrM<  JO N 


XI  101) ,Z(P1) ,U0<81)  ,PZ<81)  ,P1<81) ,P2181),P34101.81), 
UF (101*811 

Hpi.V<  101) ,F (101,81,2)  ,U(101 ,81  ,2)  ,V(10 1, 81,2)  * 

B( 101,81,2) 

IPPMT 

FC ( 101) »UC ( ICi) ,VC (1C1) ,BC (1011 
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16  Abstract 


Recently  a  new  procedure,  which  solves  the  governing  boundary-layer  equations 
with  Keller's  box  method,  has  been  developed  for  calculating  unsteady  laminar 
flows  with  flow  reversal.  In  this  report,  we  extend  this  method  to  turbulent 
boundary  layers  with  flow  reversal.  Using  the  algebraic  eddy  viscosity  formu 
lation  of  Cebeci  and  Smith,  we  consider  several  test  cases  to  investigate  the 
proposition  that  unsteady  turbulent  boundary  layers  also  remain  free  of 
singularities.  We  also  perform  turbulent  flow  calculations  by  using  the  tur¬ 
bulence  model  of  Bradshaw,  Ferriss,  and  Atwell;  we  solve  the  governing 
equations  for  both  models  by  using  the  same  numerical  scheme  and  compare  the 
predictions  with  each  other. 


The  study  reveals  that,  as  in  laminar  flows,  the  unsteady  turbulent  boundary 
layers  are  free  from  singularities,  but  there  is  a  clear  indication  of  rapid 
thickening  of  the  boundary  layer  with  increasing  flow  reversal.  The  study 
also  reveals  that  the  predictions  of  both  turbulence  models  are  the  same  for 
all  practical  purposes. 
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